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I.  INTRODUCTION 


This  report  presents  the  current  status  of  the  work  under  ARPA  order 
number  2089  ,  dealing  with  weather  modeling  with  a  view  toward  climato¬ 
logical  applications.  Emphasis  is  on  the  solution  of  the  primitive 
equations  governing  weather  by  means  of  Vector  Spherical  Harmonics.  (VSH) 
This  section  presents  a  summary  of  the  report  and  a  brief  review  of 
vector  spherical  harmonics. 


Section  II  deals  with  theoretical  developments  of  various  aspects  of  the 
technique,  notably  the  "Balance"  conditions,  advection  formulas,  kinetic 
energy,  and  an  analysis  of  sound  waves,  all  from  the  VSH  viewpoint. 

Section  III  deals  with  the  modeling  activities  proper,  including  the 
derivation  of  a  set  of  equations  suitable  for  VSH  solution  and  their 
subsequent  implementation. 

Section  IV  is  a  review  of  computational  considerations  mostly  due  to  the 
transform  method  of  handling  nonlinear  terms. 

Section  V  discusses  radial  functions,  i.e.,  the  functions  used  to  represent 
the  radial  dependence  of  the  VSH  coefficients. 

Section  VI  reviews  the  results  of  fitting  VSH  coefficients  to  published 
data,  and  Section  VII  discusses  the  current  status  of  the  programming 
system. 

Equation  numbers  are  repeated  between  sections.  All  references  to  equations 
are  within  the  same  section  unless  otherwise  indicated. 


1.1  Brief  Review  of  Vector  Spherical  Harmonics 

The  Vector  Spherical  Harmonics  are  a  set  of  orthonormal  functions  of  the 
spheriral  coordinates  ©-  (colatitude)  and  (p  (longitude) .  They  are 
triply  Indexed  functions  denoted  by 


where 


J,  L  >  O,  /m|1  |3-lU  1 

"T’M 

a  vector  field  in  spherical  coordinates  may  be  expanded  in  a  series  of  (jL 
If  the  vector  field  V  is  a  function  of  the  radius,  &  ,  and  (J)  ,  then 

vO,<s<j>)-  5ZvJt(V)T  \ 


If  the  vector  field  is  a  function  of  time  as  well,  the  coefficients  Vj. 

Z  —7—  M  L 

are  a  function  of  I C  and  T.  .  To  define  the  I  ,  we  define  firnt  the 
ordinary  scalar  spherical  harmonics  (SSH)  by 


r,  ,,  N  y/i 

,  ,  v"  CSi*lXL-*C)!  .  c-M« 

/L  (o^-c-i)  tic***)* 
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where  fl  C^icSl)  ^are  associated  Legendre  functions. 
Basis  Vectors  ^ ,  £^>y 

A  A  -A-  f  A  •  A 

e_4  =  e  *  / 


Define  a  set  of 


e,  - k 
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W^ere*  ^  are  the  usual  cartesian  base  vectors  and  i  ■  \i ~  l7 

The  vector  spherical  harmonics  are  defined  by  the  following  expression: 

_ AA  »  ^  ^ 

*>  -  Z  <  L  »■ -A  1^1  X  an  >  y" ?  ^  -  (! 


The  quantity  /A  yt*-  ^  /»■  called  a  Clebach-Cordon  coefficient. 

It  is  a  number  depending  on  L,  M,  and  Aa.  .  The  properties  of  these 
^efficients  and  of  'y^  (  o-}  <$> )  restrict  L,  J,  and  M  to  integer  values 
su-.h  that  J  >  0;  L>  0,  and  L  -  J-l,  J,  J+l;  and  |m[<  J. 

There  are  then  three  sets  oF  functions  which  we  label 
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_ rA 

The  orthonormality  properties  of  the  1  y  (__  are  expressed  by 
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Where  *  denotes  complex  conjugation  and  S  • '  is  the  Kroeneker  delta. 

J 


(9) 


The  differential  vector  operations  satisfied  by  the  I  yL  are  given 
below.  Note  that  no  derivatives  with  respect  to  &  ,  (p  appear.  Let 
f(r)  be  an  arbitrary  function  of  r  and 


=  /±H 

L~  V  3L+L 


(10) 
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Gradient : 


V-?MX/  ^\+i 


Divergence : 
.  VA 
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Product-type  nonlinearities  are  handled  by  the  well-known  formula  for 
scalar  spherical  harmonics  : 
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For  the  term  ^V,  which  occurs  in  Hydrodynamics  and  other  fields 
(V  is  a  vector  field),  the  expression  (apart  from  normalizing  coefficients) 
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For  many  applications  it  is  convenient  to  use  linear  combinations  o f  ( j 
which  are  closer  to  the  physical  coordinate  system.  Such  a  set  is  given 
by  the  set  A*  (»,<?)  ,  ^ ,  and  £  ^  Jefined  by 


,  tA  .  v  ,  - />! 

A<-  A,<p;  -  (O 
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(20) 


(21) 
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It  may  also  be  shown  that 

a7  =y1.,"<'s>,<p)  <2^ 


(23) 


(24) 
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rv\ 


(25) 


/-* 


Where 


K  -- 


(26) 


From  the  above  equations,  the  A  v)  are  along  the  radial  direction, 

while  &<_  (©,<9 )  and  are  tangential  to  the  surface  of  the 

earth.  For  the  case  M-o  »  fB>  “^fldional  and  c;c*)  is  zonal. 

In  analogy  with  equation  (2a) ,  the  vector  field  \y  can  be  expanded  as 
follows: 


The  spatial  derivative  operations  in  this  representation  are 


Gradient: 


Divergence : 

V-  P(.  r)  a”  -  Oa  fit A  V, 
V.PW6:-- 


(28) 


(29) 
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Here,  and  usually  In  the  following  equations,  we  will  omit  arguments  on  the 
functions. 

These  results  are  discussed  In  previous  reports  in  considerable  detail. 

The  advantage  of  VSH  in  solving  partial  differential  equations  is  that 
the  angular  part  of  the  equation  can  be  separated  out  using  (28)-(34); 
reducing  the  equations  to  partial  differential  equations  in  the  variables 
y  and  t  .  By  expanding  the  coefficients  in  suitable  functions  of  f  , 
the  equations  can  be  further  reduced  to  ordinary  differential  equations 
in  time.  This  is  made  possible  by  the  fact  that  none  of  the  spatial 
differential  operations  (curl,  divergence,  gradient,  etc.)  involve 
derivatives  with  respect  to  <9-  or  <$>  ,  and  by  the  orthonormality  of  the 
Vector  Spherical  Harmonics. 


Computationally,  it  is  advantageoua  to  use  the  Transform  method  to  evaluate 
nonlinear  terms  (see  report  for  previous  year).  In  this  method,  the  fields 
are  converted  from  spectral  to  physical  form  by  evaluating  equation  (2a) 
or  its  scalar  analog  for  each  variable  in  the  nonlinear  term.  The  result 
is  the  variables  represented  on  a  physical  grid.  The  nonliroar  operation 
is  then  performed  on  the  physical  grid  and  tiiC  result  converted  back  to 
coefficients  by  equation  (2b).  In  practice,  this  saves  considerable 
amounts  of  computation. 


II «1  Balance  Conditions  and  Truncated  Spectral  Forms 

The  purpose  of  the  following  sections  it,  to  investigate  the  various 
meteorological  "balance"  conditions  when  the  equations  are  cast  in 
spectral  form.  As  will  appear,  great  care  must  be  exercised  in  orde.. 
avoid  inconsistencies  or  anomalous  results. 


II. 1.1  The  Hydrostatic  Balance  Equation 
The  hydrostatic  balance  equation  Is 

t  =  -n  «  *9£  --  -  s  (l.D 

To  Illustrate  the  difficulty,  suppose  we  expand  ^  and  *  In  scalar 
spherical  harmonics  with  only  two  terms  (to  simplify  the  algebra),  then 


CY  =  -t 

«°t  Y° 

(1.2) 

¥  =  P°o  Yl  + 

J  r  .  1 

0°  Y° 

(1.3) 

the  product  is  then 

»£  .  W.r.-xipYtn-KlfilW 

d>r 

Now  by  the  product  rules 

y°0  y°0  =  y; 

Y*  Yl  =  c<  Y° 

Yf  Y"  =  c,v;  + 


(1.5) 


so  the  equation  (1.1)  is  then 


*  <pnx*cx)  d.6) 
=  ^ 


or,  collecting  coefficients  and  equating: 


(i.7) 


+■  =  o 


=  o 


(1.9) 


Now  in  order  to  satisfy  (1.8),  either  ^  or  must  be  zero  (or  both). 

If  «(*is  zero,  then  (1.7)  forces  either  <<*  to  be  zero,  e*'®  cannot 
be  zero  since  in  (l-?),  *  £  0, hence  c/^  =  0  .  This  leaves  us  with  the 
equation 


(1.10) 


If  g  is  Independent  of  r,  as  usually  approximated,  then  the  product  of 
and  ^ust  be  Independent  of  r;  that  is 


2 

c>>~ 


~  O 


(1.11) 


or,  substituting  »  ^£®  we  have  from  1.11 


<9o_ 

£  r> 


dr 

Ctr 


i  * 


^2 


-  <9 


(1.12) 


Thus  we  are  left  to  infer  that  if  g  is  Independent  of  0  and  f  ,  and 
°(  and  ^  are  represented  by  truncated  series,  then  *  and  f  are  independent 
of  9  and  .  Furthermore,  the  coefficients  of  r  and  p  are  related 

by  (1.11). 

Although  we  were  led  to  this  conclusion  by  an  example,  and  not  by  a 
formal  demonstration,  "he  difficulty  will  persist  as  we  taka  more  terms, 
l.e.,  pressure  and  density  will  still  be  Independent  of  9  and  Cf  . 


It  is  only  with  an  infinite  number  of  terms  that  (1.1)  can  be  satisfied 
at  all  (because  then  g/t  has  an  expansion).  In  fact,  of  course,  and 

are  dependent  on  Q  and  .  This  points  up  the  fact  that  constraint  (1.1) 
is  very  difficult  to  satisfy  in  spectral  form,  at  least  with  finite 
expansions.  If  g  is  allowed  to  be  a  function  of  &  and  ,  the  difficulty 
will  be  removed  up  to  the  number  of  terms  in  which  g  is  expanded. 


The  main  problem  arising  with  the  preceding  equation  ic.  the  fact  that 
we  have  not  been  consistent  regarding  the  order  of  expansion.  Thus, 
while  we  started  out  with  expansions  of  order  L=^-  ,  we  retained  terms  in 
the  product  of  order  L=  Z  .  When  dealing  with  truncated  expansions,  it 
is  essential  to  truncate  the  product  at  the  original  level.  When  we  do 
this,  equation  (1.6)  becomes 


(*>; a y;  f 


r  V* 

~i  '  o 


(1.13) 


-  Y  0 

—  3  *  T  o 


Now,  collecting  coefficients  leads  to  the  equations 


d°0  p>°,  C,  i-  =  o 


(1.14) 

(1.15) 


So  that  now,  if  (say)  are  given,  we  can  solve  the  system  (1.14)  and 
(1.15)  forp®  and  p*.  In  this  sense,  equations  such  as  (1.6)  are 
consistent,  in  the  sense  that  given  the  coefficients  of  either  o(  or 
,  we  can  solve  a  system  of  equations  for  the  coefficients  of  the 
other.  The  system  will  be  consistent  only  if  care  is  exercised  with 
the  truncation  limits  when  the  product  is  formed. 
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/<)>  lL  8  ^  (  o/l  jjL_ ,  (LO  lof  L-t  o>  .  [Y'-l-i  <L  0  1°  Itt  I  0>  .  /  I  f.,  ) 

K  'V 

(1.20) 
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<  (1.21) 


fix'"0 


From  (1.20)  and  (1.21)  we  see  that  only  the  terms  In  >'  of  the  type 
C°LC°  will  produce  terms.  In  order  to  see  the  structure  of  the 
resulting  terms,  It  Is  convenient  to  vrlte  equations  (1.20)  and  (1.21) 
In  matrix  form.  The  coefficients  will  be  denoted  by  ,  as  their 
form  Is  not  Important.  The  equations  are  carried  out  to  l  -S  for 
purposes  of  Illustration. 

I „  to  N  1°  L°  1  r 
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(1.22) 
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Since  on  the  right  of  equation  (1.19)  we  have  ^r,  we  will  have  terms 

on  the  right  of  the  form  VL  r  /  |«  ,  where  r‘*  are  the  coefficients 

of  r.  Then  upon  collecting  coefficients  of  &  ^  and  C HL  on  both  sides  of  (1.19) 

and  equating,  the  following  system  of  equations  is  obtained,  where  ^  ^ J  : 

ZrJi. 

*1CS  =  M’? 

o'lC0!  +•  c*3<-|  - 

o<,C®  +  °<SC<  =  r3 

.  rO  -i-  _/  r°  =  ^  > 

°4 C  3  ^  °'7C5 


^  C4 


-  O  (1.24) 

r  0 

t 

^  ,6  *  ? 

“l 4  ^3 

+  o'lSbj  -  0 

1 

As  we  can  see,  this  system  does  not  couple  o°L 

with  and  the  two 

can  be  determined  separately,  if  they  can  be  determined  at.  all.  Ve 
immediately  observe  some  major  difficulties.  Considering  first  the 
equations  for  the  b°L  ,  eq.  ,  (1.24.6)  forces  0 •  Then  from 

equation  (1.24.8)  s o  ,  and  so  forth.  Thus,  the  even-L  b “L  =  b®* 
are  zero.  But  the  odd-lb °L  cannot  be  determined  at  all,  as  we  have: 

°S\o  ^  +  °^ii  ~  ° 

which  is  too  few  equations.  If  we  had  truncated  at  an  even-order  L  we 
could  use  the  last  equation  to  infer  that  the  highest-order  odd-L  L[  is 
zero,  but  this  is  quite  arbitrary. 

In  similar  fashion,  for  the  C®  ,  eq.  ,  (1.24.1)  yields  C\  ,  using 
this  value  in  (1.24.3)  yields  end  so  forth.  The  result  of  truncating 
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at  1=5  yields  (1.24.5)  which  gives  a  contradictory  value  for  C%  , 

In  analogy  to  the  above  discussion  for  the  ,  the  C°v  for  odd  L 
are  not  determinable  from  the  above  in  this  example.  Extending  the 
example  will  reuove  the  problem  for  the  odd  L  ,  using  the  last  equation 
of  the  set.  But  if  we  truncate  at  odd  L  we  must  in  fact  presumably 
neglect  this  equation. 

All  these  results  point  up  the  inadequacy  of  the  geostrophic  balance 
condition  in  spectral  form.  The  reason  for  the  inadequacy  lies  in  the 
fact  that  geostrophic  balance  holds  only  at  middle  latitudes.  Spectral 
forms,  however,  are  global,  and  an  approximation  that  holds  only  at 
restricted  latitudes  must  necessarily  present  problems  in  a  global  form. 
To  further  elucidate  the  point,  we  can  transform  equations  (1.19)  to 


spherical  coordinates,  yielding: 

—  2  Si-  Co'S  9  l/cp 

=  1 

(1.25) 

nf>)  99 

2  SI  CoS  9  ^ 

=  JL  dr 

(1.26) 

rep)  dep 

Now  when  M  ■  0(r  is  independent  of  <f  and  therefore  l’.°  So  we  need 
only  consider  (1.25).  Tie  factor  (CoS 9)  causes  to  approach  infinity 
at  the  equator,  unless  also  contains  a  Col&  factor.  In  a  spectral 

representation,  contains  associated  Legendre  polynomials,  which  are 

alternately  zero  and  non-zuro  at  the  equator.  Those  that  are  non-zero 
at  the  equator  cannot  satisfy  (1.25).  Thus,  (1.25)  implies  that  only 
even-order  coefficients  can  be  determined.  This  is  unfortunate, 
since  data  analysis  shows  that  the  odd-ordered  C[_  are  the  major  part 
of  the  zonal  wind.  (The  zonal  wind  is  in  fact  given  by  C°  .) 

The  preceding  discussion  points  out  the  fact  that  the  geostrophic  balance 
is  not  a  suitable  approximation  to  make  on  any  spectral  model,  at  least 
one  utilizing  vector  spherical  harmonics. 


Next,  we  consider  the  ce^e  where  M  ^  0,  that  is,  when  the  fields  r  and 
9  depend  on  p  .  In  this  case  the  formulas  are: 

(  £.  K  Ci-<  +  Go^Hl^l  ■*  ^  )  (1.27) 


Rye 


*  c*  =  -a  (  &  cuL  1  *■  G<>  cl  c l  +  &Lt,) 


(1.28) 


where 

C-  _  =  jLMlolu~l"'> 

"5T 


’o  S 


^m  /  »  /  ^  > 


0  f.  ~  ft  L  I  ^  M  I  &  l  1  ^  ^ 

«4 

Let  us  cousider  the  case  ^mix  ”  ^  ~  3  for  purposes  of  illustration. 
We  obtain  from  the  equation  (1.19)  the  expansion 

*#■**(  l  tLSu  *C*  e?  j  =  -&r  =  H/'^u  1.29) 

In  analogy  with  the  set  of  equations  (1.23)  we  now  have  the  array  of 
coefficients  on  the  left  as  follows:  (we  have  omitted  the  6+;o;- 
factors  for  clarity) 


' 


On  Che  right  we  have  Che  terns  containing  ft* .  Upon  equating 
coefficients  of  the  6 *  and  C *  In  equation  (1.28)  we  can  write 
the  following  systems  of  equations: 


0*1  b'l  +  °^7CZ 

* 3  ^  CZ 

■■  »' 

-  '3 

o/Sl>l  *  *6  ^3  +  ^7  c2 

=  0 

*S  Q  +  c[ 

- 

*16  'y  °^n  C  2 

^  (5 

c£,l  t>\  +  U>3  CS 

r2 

—  '  2 

<*rt  +  U,s  C3 

-  0 

f  ^  *>3  +  °^I7  C\ 

^3 

\  °^iS  +  *11  c\ 

=  O 

r  <**>  ^  *  «*/  c;  +  0(22 
J  0^23  b'z  •+  *2* 

C3 

Cj  =  0 

V  ^2  5  V<l  f  C/ 

=  0 

(1.31) 


(1.32) 


(1.33) 


(1.34) 


(1.35) 


S>*  =  r. 


oi. 


ll 


c\  -  o 


(1.36) 


o?-  10 


where  the  are  coefficients  of  the  form 

<4  =  -Jii— 

and  a  similar  set  of  equations  can  be  written  for  negative 
Thus  we  see  that  for  M  i  0  the  geostrophic  balance  condition  presents 
no  formal  difficulties.  Upon  reference  to  equations  (1.31)  -  (1.36), 
we  may  observe  that  the  spectral  expansion  of  r  contains  r  \  ( col  <9  ) 
factors.  Such  factors  give  a  cosine  term  when  differentiated  with 
respect  to  &■  cancelling  the  (£etf  <9  )  on  the  left  of  (1.25).  The  other 
equation  must  be  ex.  nined  more  carefully,  which  we  will  omit,  but  no 
difficulties  arise. 

However,  this  is  not  encouraging,  since  the  case  M  ■  0  cannot  be  worked 
out,  and  these  coefficients  constitute  that  part  of  the  zonal  circulation 
which  is  independent  of  longitude.  It  is  concluded,  on  the  whole,  that 
the  geostrophic  balance  assumption  is  not  suitable  for  use  in  VSH  modeling. 


Further  Development  of  the  Advectlve  Term. 


/  ”*>  V 

The  advectlve  term (y  jV  has  been  derived  before  and  a  typical  term  fas 
the  form 


(2-1) 


l-i'-C*  i 

.  2  U’  l'  C  T  1 J  ^  1  h  M1 1  3"  fi'  )  ^-1  \  *" 

O" 


When  (2-1)  la  actually  computed  there  appear  some  simplifications  which 
do  not  seem  to  be  derivable,  at  least  not  with  a  reasonable  amount  of 
labor,  from  the  properties  of  the  various  vector  coupling  coefficient. 

In  order  to  Investigate  these  simplifications  we  have  employed  a  different 
method  which  we  will  now  describe.  The  Idea  Is  to  try  to  capitalize  on 

-r*h 

the  fact  that  In  addition  to  the  ^representation  to  which  is  orthogonal 
In  the  function  sense  one  can  use  the  representation  in  terms  of  the  linear 
combination 

chu(»,v  ) 

discussed  before,  which  In  addition  to  being  orthogonal  in  the  function 
sense  also  an  orthogonal  In  a  geometric  sense. 


c?-12 


It  xs  this  last  property  that  we  want  to  investigate  further. 


We  repeat  one  of  the  definitions  of  these  functions: 


C  -issvw,y).. .  I  y  .e  1  <5vh7 

O©  J 


tu  =  \uuu 


(2-2) 


(2-3) 


(2-4) 


(2-5) 


In  spherical  coordinators  the  del  operator  is: 


^  a-  e.  d  *  eQ  *  e  j~.  £> 

^  <  *-  —  ^  _  _v  ,  >  — 

Q  «r  <  <de  c  f*w  ©  > 


'by 


(2-6) 


The  A's  are  radial  functions  and  the  B's  and  C's  are  tangential  to  the 
surface  of  a  sphere. 


Let  us  break  up  the  velocity  field  into  two  parts.  We  will  now  use  the 
Indices  1  and  2  to  refer  to  vertical  and  tangential  components , 
respectively. 


(2-7) 


^  *  X  +  'lx 


(2-8) 

2o\'Al^t8/y) 

(2-9) 

^  t  s>,*)**  cVlV>e^M 

(2-10) 

■sj  :  <?r0.  ,  ^  s  «ee  +  *  L  4 

(2-11) 

We  can  then  write  t  ^  as : 

VVVj'lY(vV)Y^','0  -  (5^  ^ 

(2-12) 

( ^. •  ^ ^  +  (^j  ^(\v(  MM  ^ 

K 

v  (VT>V,v 


Let  us  consider  esch  of  these  eight  terms  sepsrately: 


^  L  u  <^t\w  u  0 

,  S<L  *UJ^  e  Vhlv*\ 

U\H  vV  U  v  L  1  vJ 


(2-13) 


=_L  2.1  *  Ku^.LWiiH'hVV^  >A 

ilr  Lv>  \J  \V  /'Ti\r  . »  *.  — 


We  thus  find  that  the  advectlve  tern  of  the  radial  velocity  with  Itaelf 
la  radial.  Further  we  see  that  the'Vsf  term  in  the  general  expression  Is 
In  fact  aero  for  the  conblnatlon  of  which  is  vertical,  namely 


^  s  v(v?,’0 


(2-14) 


This  allows  one  to  show  that  there  are  additional  relationships  between 
the  coupling  coefficients  which  were  unknown  before.  We  postpone  the 
investigation  of  these  coefficients.  Now  consider:  )  \J 


(i-  m  - 


L  '  L  Lv 
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<+-m  on  mi  ■ 


Thus  we  have  the  result 


Consider  ( .  \  \f 

2<<^w;(vRi 

H  *V- %\  w\ e, J  JL  j2  i  #(  «p) 


*  t  tKi<^\Ch 


u'  e^)) 


He  see,  then,  that  (M  ■''vp^is  tangential  and  again  the  /  v 


(2-15] 


-3 


(2-16 


terns  are 


missing. 


Consider 


V2  is  tangential  so  W-v  )  has  the  form 

l«\  V^t\v?^J<i8iViVyy£,  £r 

1  *  A  \ 

'4r  ;?U  0  -V^ 1 

Consider  ^ 

I 

^'V  -‘(Z  b*M  CielVl,2c^lT)c^,Y)) 

.ll'e,,?.  *  «J#  l  IfoT'n,,:  wh...  ,)) 


”  j*  <  1(22  cx*l  u\ 
*  L  ^6  >Vf  Jp  i' 


'  *  h 

e  h 

V  L, 


^If' 


^iv  *  ai * 


(2-17) 

(2-18) 


(2-20) 


Thus  this  expression  has  both  radial  and  tangential  components  but  no 
derivatives. 


_r> 


Fianlly  consider  (\/  ,  ^  ^ 

1  V 

^ ^  u"^  ^  v  -  "2  iw^a4;.  «#,y,\*,2cVl-4|  *.*, 
•  C^|v,4\  Cy%!  <P) 


»V 


(2-24) 


A- 17 


This  gives  radial  and  tangential  fields.  Since  one  can  expand  the  B'e, 

C'a  in  T'a  and  then  only  have  derivatives  with  respect  to  r  and  there  are 
no  derivatives  with  respect  r  In  the  formulation  above  it  means  that  there 

can  be  no  derivatives  at  all  in  the  general  formulation  of  this  term. 

We  see  then  that  the  only  terms  in  containing  derivatives  are 

all  the  other  have  a  factor  Vv*  ant*  r  taken  to  be  the  radius  of 
the  earth  one  sees  that  there  contributions  are  small. 


We  will  now  consider  how  the  information  obtained  above  can  be  applied 

,3  ""  '* 

to  slmplylfying  (1)  and  also  to  the  derivation  of  terms  like  W  ■s'J  V 

i*  V 


First  consider 

K  ch  t\\W  C-\|  If ))  ■  V  )  c'  £-hJ  V 

which  is  the  simpllest  of  the  terms  involved. 

Chu  £  T£u<<>. 


(2-25) 


1  u  ru" 


Since  the  d/dr  terms  vanish  one  finds 


(2-26) 


c^18 


This  result  has  also  been  derived  from  the  definitions  of  the  Z-coefficlents. 


The  'I//  terms  then  become 

LL  L  L  •  W."  )l\>ll -«/  <  uu'  +  i  <2"27) 

multlpled  by  the  terns  outside  of  the  curly  brackets  In  equation  (2-26) 

The  other  terms  can  be  derived  In  a  similar  manner.  We  will  not  display 
he  results  here  since  they  are  obtained  from  a  straightforward  manipulation 
completely  analogous  to  the  derivation  of  (2-26). 


II. 3  Kinetic  Energy  and  Angular  Momentum  in  VSH  Spectral  Forms 
The  section  on  angular  momentum  is  taken  from  a  r revious  report  (1972) . 
It  is  included  for  completeness. 


II. 3.1  Kinetic  Energy 

The  kinetic  energy  of  the  atmosphere  is  given  by  the  (.xpression 


IT  =  Jf  7.  v  cro 


(3-1) 


where  7~  I3  kinetic  energy  and  dT  is  a  volume  element.  The  Integration 
is  over  all  the  atmosphere.  As  usual,  we  expand  v  in  vector  spherical 


harmonics: 


k  »  2.  aLr  L  f  ^  '  C  L 


(3-2) 


We  will  assume  th  it  j)  is  Independent  of  0  and  (  and  a  spherical 
earth,  then 


r  „i,T  ^  IT 


(3-3) 


rr  r 

-  r  =  T  j  /  I  yV  J  ( +  tc  *  )  . 

rd  { J>  i  <<■<:+  +<:W)  *.£?.. 

The  expressions  for  the  dot  products  of  A?  and  can  be 


represented  as  follows: 


„  .V  Kt  '  S’  f1rl  .  A-f  M  , 

Al-Al'  "  IC«aYL»  (3_4)  *  0  (3'5) 


Cd  =  -  CHL-CUJ  =  l  Cei 

Sl!.c'i:  =  -^r .  J  c*c  y^' 

where  the  ,  etc.  depend  on 


M  r 


(3-6) 


(3-7) 


^  20 


In  addition,  the  complex  conjugate  relations  may  be  shown  to  be 


—  M 
*  , 

~  (-)  M  A?* 

(3-8) 

t  / 

8~i 

=  (-V  8** 

(3-9) 

=  (->*"  cl* 

(3-10) 

We  also  will  require  the  fact  that  when  /. s  L'  and  Hz  -M ,  the 
coefficients  ,  C/?n  ,  etc.,  reduce  to 


Cm*  Lit  ■■  L-)H//W 


(3-11) 


Upon  inserting  the  expressions  (3-5)  through  (3-7)  in  equation  (3-3)  we 
obtain : 


l,T  TP  / 

f  J~  f'  r  1  /  *  K  IjMl'  -  ceici  drdcfSf- 


ir-  L 

ML 


-  0  c(& 


ML'  £  0  0 


l " 


(3-12) 


The  remaining  cross  terms  will  cancel  due  to  (3-7) .  We  may  now  express 
any  term  in  (3-12)  as 

VT  T 


J  f(r )  dr  J  jf  Y  d 


(3-13) 


o  a 


r  ,  _ 

and  note  that  this  integral  is  zero  unless  M  =  -tf  and  L  -  O .  Since 
jL-L'l  <  l"  S.  L+L'  we  require  L-  Ly  and  H'z.  -M  and  all  other 
terms  vanish.  Equation  (3-12)  is  then 


0,1  IT 


2r  =  ^  ^ 


ML 


Y i  d(f 


Sin  (9  o& 


(3-14) 


0  0 


9r  21 


where 


It 

=  Jf  Clef  y'“J) 


When  we  use  (3-8)  through  (3-11)  this  reduces  simply  to 

*r  =  1  [f{r)[6W*6X’*cWJe'r  <>15) 

t 

Thus,  the  kinetic  energy  reduces  to  a  single  integral  over  r.  All 
terms  carry  energy.  This  relation  is  analogous  to  other  spectral 
expansions,  where  the  energy  is  proportional  to  the  sums  of  the  squares 
of  the  coefficients. 


II. 3.2  Angular  Momentum 

This  section  is  condensed  from  the  report  for  1972. 

The  angular  momentum  is  defined  by 

1  =  JF >7  7 

as  before,  we  assume  a  spherical  earth  with  density  spherically 
symmetric.  Expanding  v  using  (3-2)  we  have 


(3-13) 


L  - 


A  ^  n.  r  -  i.  /)*-  /  M  .  c  v  (f  ^  ']r/“T 

rcr  *  f(>)  2i  i  ^  4  J 


(3-14) 


Evaluating  the  integrals  yields 


=  xfe“A1 f *7*  <3-15> 


using  definitions  of  A,  B,  C  gives 


.  2  h<  +/">r  <3;16> 


=  /  0  + 


fft  (r )  dr  Jc'tdsi  +  j  F(,lClrf ^Lct51 


u  0-17) 

.•  “  j 


o?-  22 


-  terns  in  C*  withY®  +  terms  in 

=  r!  ,  &°f ,  &\ 

Thus,  the  only  terms  carrying  angular  momentum  are 


c 


_  I 


) 


c 


+1 

/ 


$1  with  V 


Other  terms  may  carry  angular  momentum  locally,  but  must  cancel  in  the 
integrals. 


A 


II. 4 


Analysis  of  Sound  Waves  in  VSH 
This  section  applies  Vector  Spherical  Harmonics  to  sound  waves.  No 
startling  results  are  discovered  but  some  light  is  shed  on  how  sound 
waves  behave  in  terms  of  VSH. 


H.4.1  Derivation  of  the  Sound  Wave  Equation  from  the  Hydrodynamic 
Equations 


Following  Landau  and  Lifschitz 


(1) 


,  let 


(’=  p+ft  ■  ft  ;  iW'  >  w-i: 

(4-2; 


The  continuity  equation  is  ^  =  (P  ^  ^ 

at 


or 


2  (faff1)  -  il'  =  V‘[(f0rf')f  fi  (4-3; 

^  at 


The  equation  of  motion  used  is 

+  (l /.  7;  w  =  - 

<*t  ? 


(4-4) 


or  approximately,  neglecting  (V.p)v  and  using  (4-1) 

s_j  k  ~  -Vf1 

*  P+f°  A 

We  next  eliminate  j>'  by  using  the  assumption  that  the  process  is 


(4-5) 


adiabatic.  Then 


Where  s  denotes  the  entropy.  Thus 


p'  =  (Of.  \  i;,/  or 

‘  (jfo  jS  1 

f'=  *f  '  *s  (%), 


The  equations  for  sound  are  then: 

dj_'  -  kf0 

dt- 

dv  _  L  7  p< 

*  ~  '/• 

II. 4. 2  Analysis  in  Terms  of  Vector  Spherical  Harmonics 
First  we  show  that  the  velocity  J'  has  no  C ^  terms. 


i ;  =  v<pf7x/) 

^  £  V  *  A  ( 

Obviously  the  part  of  ^  due  toV0haa  no  C^s .  For  the  -part,  let 

t ‘ia‘i  n  +  C*ci  ( 

4  =  I*  t  /)<  *-/M*  C  <■  ( 

So  (10)  is  now 


£  (  V  u-;  t  Va/1  J  -  V7  f 
Ot  Ju 

i  V<p  i-  )  -  f  V 

JQ 


(4-15) 


(4-16) 


Then  (10)  is 


ojo  'dy-i_ 


-L  7f 

f-  ! 


J  c  Jc 

Where  we  write  ^  for  convenience  (we  really  mean  f> ' ),  and  (4-9)  is 

°±  =  --  />, 

Jc  J 

Taking  the  divergence  of  (4-24)  we  have 


(4-24) 


(4-25) 


Z  (?■*,)  =  -- 

dt  * 

From  (25) ,  \7  -  j__  '  »  80  substituting  into  (4-26)  gives 

2  4 a- 

L.  .££  ~  - 

Jt  (  */u  J  if  j*J 


(4-26) 


or 


rjL  =  «  ~V 


(4-27) 


This  is  the  scalar  wave  equation  for  ^  . 

Now  we  take  the  gradient  of  (4-26) : 

2  (*f)  r  tpa  V(  7.^  ; 
if-  i 


(4-28) 


Now,  since 


7  2  l/2  =  F  (  7. 1^2 )  -  Vx7Xvt 


and 


and  since 


V  x  v~i  ~  o 


(4-29) 


V  =  U  )  from  (4 


-24) ,  (4-28)  becomes 


£  f  -+  d  V2 

dt  1  dt  jt 


k  V2? 


(4-30) 


9<  28 


Taking  the  curl  of  (24)  we  find 

9  ,  _  * 


iv  xv,)  =  o 


(4-31) 


Thus  ^i,"  const,  in  time 

or  v,  -  const,  in  time  since  curl  doesn't  affect  time 


behavior 


Chus  (30)  is  now 


<97  'vz  _  k  V 2 


(4-32) 


and  ,  the  "incompressible"  part  of  ^  does  not  support  sound  waves. 
We  now  solve  (4-27)  in  spherical  coordinates.  Let 


1  <¥ute>*>  Y:<s, 


(4-33) 


(4-33)  .ub.tltvMd  In  (4-27)  give.  equation  f<)r  the  cocfflcle„ts  . 


f  -  *  (  f,  fie-  inf.;  Ur 
Jr  \  *>  r  c)'  "r5 "  y  ' 

e)  -  <$*(') 


(4-34) 


(4-35) 


(4-34)  - 


fill")  =  2- J -*£*>]<?«  <fZ  , 

[''Jr1  r  Jr ■  f2  J‘r“  '44  v  - 


c  _  f  tl  (>•>  =  ( <?*  f  m 

ftiO)  (<"37> 


(4-38) 


29 


(4-39) 


[  V  L Hi1  + c  1  ^ 


«)  =o 


Tha  solution  to  (4-38)  is 


K,  t  T  M  -V-*c  t 

(t)  u6  *  §tiL 


(4-40) 


and  the  solution  to  (4-39)  is 


Co  = 


(rcp) 


(4-41) 


Jl  j  Vt-  are  Spherical  Bessel  functions  of  first  and  second  kind 


respectively. 


th 

We  can  also  solve  (4-32) .  It  is  more  convenient  to  use  the  fJL  rather 

than  A*  Ch  since  the  former  are  eigenfunctions  of  Vz . 

L  )  W  i  i  ^  L 


When  we  let 


we  obtain 


■sr  _  K4  -t-ki 
-  ^  ^  jL  ^ 


(4-42) 


9  V",  C,  +  9]  tr?L  r £  f  »l  C  =  Om.  'll  c 


thllh  ] 


(4-43) 


The  equations  for  the  coefficients  are; 


Jf1 


(4-44) 


c8-30 


TZ  H 

o  V_U 
'dtl 


2  ,  « 


-  -  £  f  <  4-t  4^  +  fil  An  D -i  )  V?c 


u  ft 


^-iA+i  V"UH 


(4-45) 


We  let 


C  -  ^ 


(4-46) 


and  obtain  the  following  separated  equations.  The  constant  of  separation 
18  (-1 >  Ci  >  CUi  respectively;  dr  Cj  in  general. 

We  know  in  fact  that  the  C ^  terms  will  not  be  present  from 

equation  (4-20)  (or  will  be  constant) . 

/w 


•i  ;«.<■«  +  *cjk =  ° 


(4-47) 


and  for  the  radial  parts: 

-a 


(  <L  4-  2.  d  +  L  (-Ltl)  ±  C  \ 

'  Or7  r  dr  r7  L''J 

(  il  +  li  -  i(iH) 

(  drz  r  dr  ‘  f 


H 

Al-I 


=  o 


if 


r* 


ru  O 


(  2.  +Ll  .....  !Lt:)(L*2)  .  ■  V-  4  ,  _ 

(  D'1  r  Jr - - -  -ltl  )  fiUI  - 

)r  * 


(4-48) 


(4-49) 


(4-50) 


and  therefore: 


_ .  ■  i 


Ut 


'/£>  t  ,  M 

t  +  V  tJti 


-  \J-fcj  t 

£  (4-51) 


^31 


irAj  —  \J 

w  '  ■  -  iruL-t 


(4-52) 


'V 


i,  +'C-,  yt-/ 

r  C,  jL  (V^"f)  ■+  '/2Ml  j/t  (  )  (4-53) 

*  /VI  ^  ^  '  f 

r-‘-"  =  ^  zr^,  [ju,  i  /»)  <4‘54> 


/■  .-/ 
N 

>‘L-L 


«y  At 


Where  CL  Q^,  must  be  determined  from  boundary  conditions;  and 
are  spherical  Bessel  functions  of  the  first  and  second  kind.  Since  in 
fact  we  know  that  TfL are  not  present  we  can  dispose  of  (4-53)  in 
considering  wave  motion. 

11.4.4  The  Velocity  Potential 

Since  the  1 LL  j  are  not  in  fact  present,  we  can  write  v  as  the  gradient 
of  a  scalar  velocity  potential,  paralleling  conventional  developments. 
(Note  that  it  isn't  necessary.)  When  we  do  this  we  have 


v  =  V  ^ 

/•> 


3J  =  -:< 

dt 


(4-55) 

(4-9) 


Then 


<3_2  i  'L  Vf  (4-10) 

a>t  fa 

tC  =  -k]o  (V-  (V*p))  -'kJ' 


dt 


d(V\b)  _  ¥(  djt  )  =  -  L  V*p 

7T  ’  {'tJ 

So  that  aside  from  a  proportionality  factor, 


t 


(4-56) 

(4-57) 


(4-58) 


cf.  Ref,  (1) ,eq.  63.6.  Using  (4-58)  and  (4-56)  we  get 

dj  -  =kPt,VZf 

*  t  i)  c  7- 


or 


32 


Which  is  identical  to  the  pressure  equation  (4-27)  so  its  solution 
is  given  by  (4-40)  and  (4-41),  with  different  constants. 


II. 4. 5  Boundary  Conditions 

Note  that  in  (52)  and  (54)  we  have  the  constants  CL_f  f  Cltkt.  We  have  to 
infer  those  constants  from  boundary  conditions.  If  the  boundary  condition 
is 

v.  4  =  o  ,  =  p0 


that  implies 

aHL(*0) 


=  ft 


Un-t  ~  ^  U 


LLt! 


V, 


V  =  O 


(4-60) 


or 


Pl  V 1KL-I  Jl-I  ro)  Vj rLu,  jLfl  A,  )  ~t 

A  ^Au-/  -a-  =0 


(4-61) 


This  yields  one  equation  relating  .  With  the  boundary  condition 

chosen,  cannot  be  determined.  This  conforms  to  the  fact  that 

the  wavelength  (or  frequency)  of  sound  waves  on  the  surface  of  a  sphere 
have  no  characteristic  frequencies.  If  we  were  to  impose  an  additional 
condition,  e.g.,  V  •  tr  r  O  when  I1:)]  ,  we  would  have  an  additional 


equation,  similar  to  (4-61) ,  relating 
the  characteristic  frequencies. 


O  ,  C  ;  which  would  determine 
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-II*  Prioitive  Equations  Meteorological  Modeling 


T':1S  “Ctl0n  dl8Cuase“  eff«“  develop  .  numerical  weather  prediction 
model  based  on  a  dry-air  verelon  of  the  Mlntz-Arakaw.  Model  «>  uelllll 

the  method  of  Vector  Spherical  Harmonic,  for  the  solution  of  equation.. 

*  previous  report'2)  described  a  quasi-an.!,^  solutlon  „f  t„e 
For  the  current  effort,  the  level  of  Implementation  Is  confined  to  the 
framework  of  the  General  Program  System  described  In  the  same  report.  The 
modeling  activity  serves  as  a  test  for  the  general  system. 

In  analyzing  the  Mlntz-Arakaw.  model.  It  appears  that  It  Is  distinguished 
from  other  models  by  (1)  The  forcing  fields  assumed  by  the  model,  l.e. 
humidity  Input,  radiative  heating,  and  friction  loss  terms:  (2)  The 
horizontal  differencing  technique;  and  (3)  the  method  of  handling  convective 
adjustment,  of  these  aspects  (2)  1.  not  applicable  to  Vector  Spherical 
Harmonics  rad  (3)  does  not  apply  l„  view  0f  the  dry-air  assumption  which 
.1.0  remove,  humidity  Inputs.  Since  so  much  of  Mlntz-Arakaw.  proper  doe, 
not  apply  to  this  effort.  It  was  deemed  constructive  -  especially  given  the 
fact  that  a  novel  solution  technique  was  to  be  used  -  to  begin  with  the 
basic  primitive  equation.  In  spherical  coordinates  and  adapt  them  for  VSH 
solution.  A.  a  result  of  several  trials,  a  model,  or  set  of  equations 
almost  identical  to  Mlntz-Arakaw.  was  evolved.  The  following  sections 
describe  this  evolution  In  some  detail.  As  an  interesting  result,  the 
primitive  equation,  have  been  rigorously  converted  to  spherical  coordinates 
using  successively  the  radius,  the  pressure,  and  the  normalized  pressure  (r) 
as  a  vertical  coordinate.  In  the  process  of  conversion,  a  number  of  terms 
are  chained  which  are  usually  neglected  In  the  literature.  In  the  final 
verszon  of  the  model  described  herein,  thee,  additional  term,  have  dropped 
to  maintain  consistency  with  the  Mintz-Arakava  model. 


Section  III.l  states  the  assumptions  used  to  derive  the  equations; 
Section  111.2  presents  the  derivation  using  p  as  a  vertical  coordinate 
without  the  hydrostatic  balance  condition.  Section  111.3  discusses 
boundary  conditions.  Section  111.4  the  limitations  of  the  model.  In 
Section  III. 5,  the  hydrostatic  balance  condition  is  imposed  and  the 
changes  to  the  model  are  described.  Section  III. 6  covers  the  change  to 
pressure  and  sigma  coordinates.  Section  III. 7  presents  a  functional  flow 
of  the  implementation . 


III.l  Assumptions 


The  assumptions  used  to  derive  all  equations  in  the  following  sections  are 
o  Dry  Air.  No  humidity  terms  are  considered  in  any  equation, 
o  No  orography,  (in  radial  coordinates) 

o  Spherical  earth. 

o  The  earth  is  assumed  to  be  "all  oceans".  That  is,  the  surface 
temperature  will  he  taken  as  a  constant  function  of  9 ,  cj)  . 

In  deriving  the  first  set  of  equations  -  which  form  the  basis  for  much  of 
the  rest  of  the  work  -  the  Hydrostatic  Balance  approximation  will  not  he 

used.  Also,  in  the  initial  equations,  the  earth's  gravitational  field  is 
arbitrary. 


The  terminology  in  the  following  section  is  given  below: 


v  three-diminsional  velocity  vector  field 
f  density  field 

_  A. 

earth's  angular  velocity  vector 
pressure  field 
T  temperature  field 
R  universal  gas  constant 

K  Ratio  of  Cp/Cv 

Yc-  gravitational  constant  CM  ,  where  M 

e  e 

■  mass  of  earth;  G  ■  universal  constant 
(gravitational) 

^  frictional  loss  term  vector  field 

Me'  heat  input  function/temperature 

/V  heat  input 

(p  earth  gravitational  potential 


(meter  sec  1) 

(gm  meter  ) 
(sec-1) 

—1 

(gram  meter  sec 
(  °K) 

+2  —2 

(meter  sec  (  K) 
(specific  heats) 

,  3  -2. 

(m  sec  ) 

_2 

(m  sec  ) 

(■ec_1) 

(ly/sec) 


Unless  otherwise  stated,  all  fields  are  functions  of  V'  ,  <9  ,  <p  and 
all  vector  fields  are  three  dimensional. 


The  above  set  of  equations  is  not  suitable  for  solution  by  means  of  the 
General  Program  System;  for  which  it  is  preferable  to  have  an  equation  for 
the  time  derivative  of  each  variable  (vector  or  scalar) .  One  way  to 
obtain  such  a  system  is  to  eliminate  p  from  the  equations  by  means  of  the 
gas  law  (2-3).  Using  this  equation  in  2-1  produces: 


ft1  -  -(V-  Vp.  2J<L  X  /  p  t  F 

p 

Equation  (2.2)  may  also  be  written  (see  2-14  }as 

dp  -  —  P  V *  v 
oc 

And  using  the  gas  law  to  eliminate  dp/dt  in  (2-4)  results  in 

—  ‘■-if  -  ci~  _  krt  ctp  -  C 

T  V  <7?  7  JP  “ 

Substituting  (2-8)  into  (2-9)  and  collecting  terms  gives 


(2-7) 


(2-8) 


(2-9) 


il  -  2X  -LL  V.u  ,  ,  (cj_  tn.v) 

cit  ( >  -  /< ;  /-a'  /-a 

If  we  now  use  the  gas  law  to  eliminate  dT/dt  in  (2-4)  we  obtain 


(2-10) 


— -  '•/'  t  X  '■!£  -  Z—  dJL  -  Q  (2-u) 

p«r  'It  f2cr  Jc  f  .it 

again  using  (2-8)  in  (2-11)  and  collecting  terms  in  dp/dt,  we  find 


J  cl£  _  (0--  V'  *  ) 

t  ciC  !-K 


(2-12) 


&  6 


N 


t:  aprtral  ”eth0dS-  “  ls  to  avoid  divisions  by  s 

variable  whenever  possible.  Such  a  division  oooura  In  the  equation  of 

notion  where  Vf.  i,  divided  by  f  .  To  avoid  this,  introduce  the 
variable  ,  defined  by 


</s  2  f 


Thus  Vf/f,  =  V<^  and  t/  Jf/Jf  =  d^/dtr 

Finally  we  use  the  fact  that  if  f  is  a  scalar  field, 

>jf  =  sf 

lit 

The  complete  set  of  equations  is  then; 


°I  y-  l/.  yf 
dc 


—  -  -(V'V )  v  -  /<7'V<p  -  ?SZxv  /  ty(f>  /- 


(2-13) 


(2-14) 


(2-15) 


dt 


Jr 


~  —  *  V  tj;  -y- 


/—A; 


—  r-l/'7T  -f-  _l_ 

!-K 


(<f-  V-  V  ) 


9 

(Q  -  av-  v  ) 


(2-16) 


(2-17) 


This  is  a  set  of  1  vector  and ^2  scalar  equations  in  vector  unknowns 
f  an  scalar  unknowns  ^and  /  .  T*ese  equations  are  suitable  for  solutic 
ns  ng  the  General  Program  System.  It  is  interesting  to  note  that  (2-16) 

‘  “7"  "CePt  ,<>r  th*  '™  <»— **  small)  ?.  ,  and  that 

variabl  '  (t°  —  ■""—•«-)  by  defining  a 

variable  equal  to  („  T  .  Thi.  change,  however,  would  explicate  the  equat 


v*-7 


of  motion.  To  complete  the  specification  of  the  problem,  we  must  state 

•  * 

the  boundary  conditions,  the  form  of  the  heating  term  Q  or  /-/  ,  and  the 
form  of  the  friction  loss  term. 

To  Implement  a  model,  we  must  In  addition  specify  the  set  of  radial  functions 
to  be  used  (3ee  Section  V),  the  number  of  layers  to  be  used,  and  the  layer 
spacing. 
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III. 3  Boundary  Conditions 


The  boundary  conditions  to  be  specified  are: 

1)  The  normal  velocity  component  -  in  the  radial  direction  for  a 
spherical  earth  -  is  zero  at  the  earth's  surface,  or 

{V'?rK  =  r0  =  0  «-» 

2)  The  surface  temperature  is  treated  as  constant, 

T  (r0l  tf,  $  )  -  /~  (  &  y 

—  (3-2) 

Where  is  an  input  to  the  program. 

3)  No  boundary  conditions  on  are  specified. 

The  fact  that  a  frictional  term  is  present  really  requires  that  the  velocity 
on  the  earth's  surface  be  zero.  In  reality,  this  leads  to  boundary  layer 
considerations  which  are  best  avoided  for  the  moment. 

In  the  Mintz-Arakawa  model,  the  equivalent  of  the  radial  velocity  (  O’  ) 
is  set  to  zero  at  the  top  of  the  atmosphere.  For  equations  (2-15)  -  (2-17) 
this  equivalent  step  was  not  performed. 
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III. 4  Model  Limitations 


Equations  (2-15)  through  (2-17)  constitute  a  complete  set  of  equations, 
without  any  balance  conditions  of  any  sort.  Thus  they  are  subject  to  wave 
solutions,  that  is,  the  same  sort  of  computational  instabilities  that  give 
rise  to  the  well-known  Courant-Friederichs-Lewy  (CFL)  condition  for  the 
numerical  integration  of  hyperbolic  partial  differential  equations .  We 
may  identify  the  part  of  the  equations  that  cause  trouble  by  considering 


(4-1) 


only  the  equations 

as  follows: 

dv 

=  -*'7 

Vp' 

at 

dp 

V* 

at 

/-* 

(4-2) 


That  is,  the  nonlinear  terms  have  been  dropped,  the  temperature  has  been 
held  constant,  and  the  corlolls  force  and  all  forcing  terms  have  been 
neglected.  In  this  simplified  model,  taking  the  gradient  of  (4-2)  gives 


vf  2  (W  =  -jl 

1  c>ty  /-*. 

taking  of  (4-1)  gives 


(  )z  -  I  y 


7- 

V 


=  -kt  §  ( v<i>) 
°  at 


combining  (4-3)  and  (4-4)  gives 


(4-3) 


(4-4) 


dt2 


!-K 


v  V 


(4-5) 
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This  is  the  second-order  wave  equation  identical  to  that  studied  in 
Section  II. A  when  dealing  with  sound  waves.  In  fact,  the  propagation 
speed  of  the  wave  is  given  by 


c  = 


(4-6) 


Which  is  approximately  the  propagation  speed  of  sound  waves.  Thus  the 
solution  contains  sound  waves. 


We  can  also  expect,  from  the  similarity  of  equations  (2-16)  and  (2-17) that 
a  similar  result  holds  for  waves  generated  by  the  equation  for  /  and  the 
equation  for  T.  The  fact  that  the  equation  contains  a  T  multiplier  makes 
it  more  difficult  to  analyze. 


With  respect  to  solution  by  spectral  methods,  the  problem  is  reduced  to 
solving  (numerically)  a  set  of  ordinary  differential  equations  in  time. 
Once  the  radial  functions  are  chosen,  the  partial  differential  equations 
in  the  coefficients  become  ordinary.  For  each  coefficient  there  are  K 
ordinary  differential  equations,  where  K  is  the  number  of  radial  functions 
chosen.  The  numerical  stability  of  a  solution  procedure  for  such  a  set  of 
equations  may  be  analyzed  by  linearizing  the  system  (as  in  (4-1)  and  (4-2) 
and  examining  the  eigenvalues  of  the  resulting  coefficient  matrix.  The 
stability  of  the  technique  is  controlled  by  the  largest  eigenvalue  in  the 
system.  This  is  equivalent  to  the  CFL  condition. 


The  sound  wave  coupling  between  these  equations  is  indicated  by  diagrammatic 
representation  below:  _  _  _ 

_ _ -  ( -fiHic  / o* )  -  ~  \ 

_  _(j7, 1 

dc  :  ' 


_ (Pne  inreot  ^T/a^y 


.  .i 


r 


_  _  7.7^  +  -L-  (i -[?***) 

- - -  —  /-/t 

dt 

li 


_ 


with  a  sirailiar  coupling  for  the  T  equation.  This  type  of  coupling  always 
leads  to  oscillation,  since  it  in  essence  corresponds  to  a  simple  harmonic 
oscillator  (or  a  set  thereof) . 

In  order  to  break  this  coupling,  it  is  customary  in  meteorology  to  resort 
to  the  hydrostatic  balance  approximation.  The  problem,  of  course,  with  any 
sort  of  constraint  upon  the  equations  is  that  information  (possibly  meaning¬ 
ful  solutions)  will  be  lost.  It  should  be  noted  that  other  possibilities 
exist.  For  example,  in  the  second  equation  above  one  could  neglect  the 
term  F*  v  .  This  is  not  equivalent  to  setting  7 (incompressible 
fluid)  :  but  rather  in  assuming  that  the  term  is  negligible  in  the  equation 
above.  In  terms  of  vector  spherical  harmonics,  this  means  that 


Since  (\L  is  very  small  'Q  (  (O'7  -v.r  <xT  ')]  and  /'  very  large  ,  Q  I/O  G  M  \ 
f.hen  this  implies  that  the  assumption  is  that  *  y.  ^  da?  .  0 

because  d L  are  0(X  *isec~'),  that  is  that  the  vertical  derivative  of 
the  vertical  velocity  is  small.  This  assumption  must  ue  tested  experimentally; 
but  the  above  discussion  makes  it  appear  at  least  plausible. 

Other  means  of  avoiding  the  coupling  above  are  under  investigation.  It  is 
perhaps  possible  to  devise  a  suitable  integration  scheme  or  to  find  better 
approximations. 
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III. 5  The  Hydrostatic  Balance  Condition 


The  numerical  instability  of  the  complete  set  of  equations  (2-15)  through 
(2-17)  have  led  meteorologists  to  impose  the  Hydrostatic  Balance 
Condition,  which  in  cartesian  coordinates  is 


%  -  -n 

or  in  spherical  coordinates, 


In  the  equivalent  Set  of  equations  (2-15)  -  (2-17)  this  is  also 
equivalent  to  the  equation 


*r  =  (V*  )t  (5_3) 

dr  0 

Where  the  subscript  f  denotes  the  radial  part  (along  )  of  , 

which  is  usually  approximated  by  the  constant  acceleration  of  gravity  g. 


III. 5.1  Effects 

The  first  effect  of  the  hydrostatic  balance  condition  is  to  break  the 
coupling  between  equations  (  2-15  )  and  (  2-16  ).  This  is  because  the 

velocity  field  is  broken  up  into  vertical  (  V,  )  and  horizontal  (  )  parts 

and  the  vertical  part  of  equation  (  2-15  )  is  "replaced"  by  the  hydrostatic 

balance  condition.  This  is  a  rather  glib  statement,  in  as  much  as  (5-2) 
is  not  an  equation  of  motion. 


It  is  not  immediately  clear  how  (5-2)  may  replace  the  vertical  motion 
part  of  (  5-3  ).  One  possibility  is  in  fact  the  set  of  equations 
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^6"  z.  '■> 


(5-4) 


cPi/i 

£>(- 


(  2  iZ.  x  1/  J  ■f*  Ay 


[li/-V)V  ] 


(5-5) 


3  (p 

Jt 


V  ‘  V^P  -A-  -i—  f  Q  -  7  •  1/  ) 
/-< 


(5-6) 


^  _  -  v-vr-t-  r_  (a-,c^,v  ) 

dt  y 


I- 1C 


(5-7) 


1/  = 


-  y/  ^  +  «2 


(5-8) 


That  is,  the  hydrostatic  balance  condition  is  used  to  simplify  the  equctior 
of  motion.  The  vertical  part  is  Integrated  explicitly,  and  the  vertical 
and  horizontal  parts  added  together  to  form  a  total  V" .  This  approach  has 
apparently  not  met  with  success,  as  it  is  nowhere  reported  in  the  literature 
thus  far  examined.  It  appears  ,  however,  that  this  approach  is  not  equivalent 
to  the  approach  taken  in  the  Mintz-Arakawa  and  similar  models,  where  the 
Hydrostatic  equation  is  retained  as  a  basic  equation.  Instead  of  an  auxiliary 
equation  used  to  simplify  the  vertical  motion  equation. 
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Therefore,  the  approach  that  must  be  taken  ia  to  retain  equation  (5.4) 
but  to  compute  a  total  V  in  such  a  way  that  the  hydrostatic  equation  (5-3) 
is  satisfied.  In  equation  (J-6)  and  (5-7),  the  term  V-  V  involves  both 
V<  and  ^  »  as  lt:  18  total  1/1  that  must  enter  into  the  equations.  Now, 

(5-3)  relates  the  independent  variables  T  in  (5-6)  and  (5-7)  Since 

VL  is  given  by  (5-4),  we  then  determine  to  satisfy  (5-3).  Schematically 


<t  ~  ^  *7,1?) 

T  ~  r  (  r,  j 

g  (^>T)  --  o 


(Equation  of  motion) 
^  -  equation  (5-6) 
T-  equation  (5-7) 
Hydrostatic  equation 


What  we  must  then  do  is  to  determine  VJ  so  that  the  last  equation  above 
is  satisfied.  We  will  perform  the  calculations  in  spherical  coordinates. 
This  will  introduce  some  differences  with  the  appearance  of  the  corres¬ 
ponding  formulas  in  cartesian  coordinates.  Wa  will  first  require  the 
expression  for  the  V  operator  in  spherical  coordinates. 


V  £ 

er  &  t 

<L  -f- 

£ 

whence  we  define 

r  dr 

r  db 

v,  = 

and 

C  d 

Or 

1 

^  * 

l  $  9 
t  d& 

+  !  -  4  & 

ro/'ic 9 

(5-9) 


(5-10) 


(5-11) 
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It  is  also  convenient  to  regard  this  operator  as  the  result  of  a  gradient 
operation  along  a  surface  of  constant  P  .  in  this  case  Vz  will  also 
he  written  as  .  Also,  we  have 


’  ($■**)'*  + 


=  02vr  +■  V72»  ^ 


(5-12) 


Now,  using  (5-9)  and  (5-12)  in  equation  (5-6)  we  have 

—  =  -(  \e,  +■  *4  ;  •  (  9yi  +  Vjy  )  t 
Jt  dr  /_^ 


*  jl.  (4 ...  v2  ] 


(5-13) 


We  must  now  differentiate  the  entire  equation  above  with  respect  to  r 
in  order  to  apply  (5-3).  This  gives: 


-  dJt 
d,'\  ; 


)  =  ^  (dJA  =  i  f  -1  ) 

ot  1  Jr  J  dt  \  fiT  J 


_  dr 
,?rz  Jt 


9__ 

nr  dr 


9^,  _  9_  dr 


-  f  v<  -  -L\  ~^/  +1^1  (5-u> 

r  j-tc  l  dr2  Y'x  t  dr  J 


+  i 

/  -  k  dr 


■a-*)  V  Vi<f  - 
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..  ..  *  ******* 


Applying  (5-3)  allows  us  to  equate  the  right  side  of  (5-14)  with  9/rf  times 
the  rig>-t  side  of  (5-7).  We  can  also  apply  (5-9)  and  (5-12)  to  equation 
(5-7),  We  obtain: 


-  ±  f-  fj (j±  >  - 1  7  JL 

/-<  J^i  r  J  !-*■  dr 


T-  2-  dI o-fc)+  \  ft  +  —i  [(*-, 

-  dr  L  1-k.drL 


')  K- 


—  F-v,  Jr 

erl  L  ^ 


trr 

(/-<) 


[(ic-D  r  _  v2.  vz  -t2  __ 

(5-15) 


We  may  now  collect  the  terms  in  ^  to  obtain  an  equation  for  ^  which 
will  cause  the  hydrostatic  equation  to  be  satisfied;  namely: 

ay,  +  (1-1)  +  (-Z-  t  *iK 

at*  K  *  «t)  dr  {  rz  er) 


■  -  f (<r<)  »r  V2^  -  Vz.V2+4  J  t  [(t-nfy \T-  V7'?z  i-  Of  J 

dh  *• 


(5-16) 

This  is  a  very  complicated  second-order  differential  equation  for  Vt 
Although  linear  in  l/,  ,  it  involves  very  cumbersome  computations  in 
order  to  find  either  the  coefficients  or  the  driving  function  (right 
side  of  (5-16)),  The  set  of  equations  to  be  solved  is  now  (5-4),  (5-6), 
(5-7),  (5-8)  and  (5-16)*Because  of  the  complications  of  (5-16)  which 
would  have  to  be  solved  at  each  grid  point,  we  will  abandon  this  set  of 
equations  in  favor  of  a  change  to  pressure  coordinates  and  coordinates 
The  appearance  of  equation  (5-16)  is  the  deciding  factor  in  the  decision 
to  change  coordinates  at  this  stage  of  program  development.  It  should  be 
noted,  however,  for  future  reference  that  if  T  is  given,  the  homogeneous 
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portion  of  (5-16)  might  be  solved  analytically,  which  in  some  future  system 
might  be  used  to  advantage.  The  computations  for  the  forcing  function 
would  still  be  with  us;  and  the  computations  would  still  have  to  be 
performed  at  each  grid  point  and  then  transformed  back  to  coefficients. 
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s**v: 


III. 6  Change  of  Coordinates 


This  subsection  concerns  the  form  of  the  equations  when  the  vertical  co¬ 
ordinate  r  is  replaced,  first  by  the  pressure  and  then  by  the  dimensionless 
variable  <r  .  In  this  section,  the  meaning  of  the  variables  IT  and  'ft  differs 
within  subsections.  However,  the  meaning  is  clearly  indicated  within  each 
subsection.  Section  III. 6.1  recapitulates  the  equations  derived  in  section 
III. 5.  Section  III. 6. 2  deals  with  the  implications  of  changing  radial 
variable  in  spherical  coordinates,  which  introduces  some  differences  from 
the  usual  procedure  on  cartesian  coordinates.  Section  III. 6. 3  deals  with 
the  pressure  coordinate  transformation;  the  independent  variables  being  p, 

6,  cf>  t  t.  Sections  III. 6. 4  deals  with  Sigma  coordinates  where  the  Independent 
variables  are  (T  ,  ©,  (p ,  and  t.  The  final  form  of  the  model  is  given  in 
Section  III. 6. 5;  and  section  III. 6.6  concerns  Model  Implementation  considera¬ 
tions  . 

III. 6.1  The  (r,  G,  0,  t)  System 

These  equations  are  essentially  a  summary  of  the  development  given  in  III. 5. 
Some  detail  has  been  added. 

Define: 


'Y'ty  \  ~  (log  pressure) 


irCe^t')  =  >Cire,©,  <*,t) 


(log  surface  pressure) 


-^-19 


9 


The  Model 
W. V  . 

it 


^TT  _ 
it  " 


IT 

it 


Is  then: 


'v.  5V1  -(v.-v^y,.  - 
■hv,  sv.e^  +K 


Z-ft-  co*  (eh*  v,^ 
(6-1) 


-  V^T  + 


(6-2) 


(6-4) 
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with  the  supplemental  Equations 


fi.  X 
r  '  kt 


the  Initial  Conditions  required  are 


Vv  Ci-,  e ,  cf,  tc) 

■v  (e,(f ,1-.) 

T  ( ►,  e,  <?X) 

and  the  Boundary  Conditions  are 

Tdrt  ,e,  tf  jt1)  =  TjCe,  cf)  = 

V-Oe)e,  <f,t)  =  ir(©,  <f,t) 

V,  O'., 6,^,1)  -  V,(«(f5<Cf,-t)  =  o 


II I. 6. 2  Implications  of  Changing  Radial  Variable 


The  "horizontal  del  operator"  (  Vy.  )  is  defined  as 


r  r 


A 

e. 


_ 


(6-8) 


for  gradient  and  divergence.  The  notation  ^  x  V 
"vertical  curl"  of  v"  . 


represents 


That  is: 

V,kV  s  {l**)«  V }]  eh  t  Ot-)*V 


Using  this  notation  the  total  derivative  operators  on  a  scalar  field  £  and 
a  tangential  velocity  field  \j  are 


cUl 

At 


V 

'  Sr 


(6-9) 


ciV% 

At 


St 


yyv 

Sr 


(6-10) 
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a 


where 


V,  -  J-t  and 

(v^-Vr-)  Vx  =  -  V,  «  (7f 

To  change  the  radial  variable  from  r  to  say  T~ 
del  operator”,  [understood  to  be  applied  on 


where 

r  -  <rO,e, 

V)t) 

]  is  defined  as 

7a- 

1. 

V 

r 

T 

5 

for  gradient  and  divergence  and 


The  equations  relating  V h  to  Vr  are; 

vTe  =|v,e  +  (?Tr) 


,  then  the  "horizontal 
surfaces  of  constant  <T, 


(6-11) 


(6-12) 
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W,  -  t  VVx  + 

(6-13) 

(6-14) 

where  C  is  a  scalar  field  that  can  be  functionally  represented  as 

e'C^O.cf^)  or  e"6r,e,«r,t)  .  Similarly  Vv  is  a  tangential 

vector  field. 

The  total  derivative  operator  becomes 

da-  O'  -7  _ 

at  ~  +  To-  31  +  T  v*'  7<rfe 

(6-15) 

d  ^  3Vj  3t  q- / —  v  — 

at  '  si  ’  s>  at  +  UV  VK 

(6-16) 

where 

III. 6. 3  The  (p,  0,  0,  t)  System 


The  equations  listed  in  III. 6.1  were  converted  to  pressure  coordinates  as  a 
step  in  obtaining  the  model  in  T  -  coordinates.  The  resulting  equations  do 
not  seem  to  be  computationally  useful  since  there  is  difficu.1  in  specifying 
boundary  conditions. 

In  the  conversion  the  more  general  equation 


is  used  in  place  of  eqn  (6-2) . 


Define 


The  model  is  then 


(vv  Vp)  Vv 


-Isic^eCe.^)  V, 


(6-18) 
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RT 


(6-23) 


/° 
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III. 6.4  The  (TJ&/  tp,  "£  )  System 

The  pressure  coordinate  equations  of  III. 6. 3  were  converted  to  {} — 
coordinate  form.  These  equations  are  listed  below. 

Defining: 


nr  == 


(6-24) 


where  is  the  surface  pressure 

and  Pj-  is  some  constant  reference  pressure 


(6-25) 
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\r 


ttRT 

If 


(6-30) 


^Co 
S  T 


^<r 


With  the 

P* 


13  \7  r?  l-  .  't'tt  RT  —  „ 
e  -•  <r  ^  <jp  V-‘  VT 

3LT  tt  .  \V  __.  Zv/ 

f  7<r  v‘+  (rf  ,v»7r"  r  V| 


<r  rr 


T  f-l  V/ 


T  -  -  *vT. 


—  V, 

l~  1 


60 

Tr<r 


•  !£ilTv  -  0-A}RTcoir 

y,  - 


%fv 


+  1^.  V  ^  - 

^  V<r 


RXtt 

*p 


Bizr  v  y 

%r  1  " 
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supplemental  equations 


it  r  -t-  pT 
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(6-31) 


(6-32) 


(6-33) 


(6-34) 


■»  v 


where 


9s  =  gO-fr-0 


^  Co 

77 


-rvjr 


V,  +  £T  Lv 

*■  I _  <T 


3<r 


VfrTT 


(6-40) 


Vt-  _  x 

dT  <T 


T.  vr*  v. 


CO 

irg- 


(6-41) 


T**  /  _ 

The  term  -^  (  Vx  •  V*  )  Vt  does  not  appear  to  be  in  the  RAND  equations, 

but  is  clearly  included  in  the  equations  presented  by  Langlois  ^  Kwok, 

Ref  (3)  .  We  are  assuming  that  we  have  misinterpreted  the  RAND 

notation. 

Equations  (6-27)  thru  (6-33)  have  been  developed  entirely  in  spherical 
coordinates. 


Note  that  the  following  approximations  applied  to  these  equations  will 
generate  the  equations  (6-  36)  thru  (6-  43)  . 


h 


„  '  yj*  MbSUMHM 
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O 


except  in  horizontal  momentum 


equation  (6-  27) 


where  V,  O 


Equation  (6-  33)  is  not  needed  since  all  occurrences  of  V,  in  the  other 
equations  have  been  eliminated. 


III. 6. 5  Final  Form  of  Model 

To  obtain  the  form  of  equations  (6-34)  to  (6-41)  most  suitable  for  a  spectral 

solution,  use  equation  (6-37)  to  eliminate  the  variable  from  the  problem. 

Then  eqn  (6—41)  becomes 


v.-vr 


j_  Mr 

it 


(6-42) 


Integrate  (6-42)  0Ver  <T  from  1  to  0  using  the  boundary  condition  <f  =  O 
at  T=  1  <r-o 

TT  Tfc  H  V  V.  ^V4.  vrir)  dr 


(6-43) 


Then  define 


'f  -  -C^TT 

<P  =  ^Ci"  !-«.) 


and  let 


pr  -  ° 


and  let  V<p  include  jL 


as  a  factor  to  obtain  the  model 


+  vT?S  +  Rr(^)vry-  = 

f  >(?)  +  =  A 

ir  *■  Cpir 

~r  (?)  +  A(?)  [H  +  T  +  *»•  vr>^j 


(6-44) 


(6-45) 


£  =  -*«) 


(6-46) 
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k,K|w  + 


/  1/  w 

ST- 


o 


in  layer  1 
in  layer  2 
in  higher  layers 


=  k4  w2(© -Ks')  m**  [w.(n.t  -  <e  +  K}  }  o]  e 


+  k<,Ts  e 


^K10e 


The  required  initial  conditions  are: 

(  <r,  6,  <p,  to) 

T  ( r,  e  , 

Y  (Si 

The  boundary  conditions  are: 

e,  ‘P.'t)  =  o 

<f  (l,  e,  <p,t)  -  r (°A  <tvt)  - 

T (l,  d.yft)  -  V®,  <p) 


O 


III. 6. 6  Considerations  for  Implementing  the  Model 


The  radial  expansion  functions  to  be  used  in  this  implementation  are  Legendre 
polynomials.  They  are  orthogonal  over  the  interval  (-1,  1).  The  variable  S 
is  a  mapping  of  <T"  that  is  used  as  the  argument  of  the  Legendre  polynomials 
in  the  implementation.  Computations  are  performed  at  five  values  of  C" 
(computational  layers).  Figure  6-1  illustrates  the  relationships  between  <T  , 
S,  p  nominal,  nominal  at  the  five  computational  layers. 

Figure  (6-1)  Radial  Structure 
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The  values  are  the  roots  of  the  equation 

P £■  CS)  -  o 

and  were  chosen  because  they  are  the  sample  points  for  Gaussian  integration. 
The  actual  mapping  relationship  is  given  by 


s  =  1  +  (s,  - 1)  <r 


(6-54) 


For  the  implementation  of  the  "horizontal"  del  operator,  where  \7j-  Includes 
the  factor  ~  ,  we  have 


<l^  0  op 


Then  the  "horizontal"  gradient,  divergence  and  the  vertical  curl  of  scalar 
and  vector  spherical  harmonics  are  given  by: 


iwW.,  (6_: 
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IJI.7  Implementation  Functional  Flow 


The  solution  of  equations  (6-44)  to  (6-48)  involves  integrations  over  time 
and  s^-gma.  The  equations  are  all  formed  in  physical  space  and  all  inte¬ 
grations  are  performed  in  spectral  space .  The  sequence  of  integrations 
needed  to  solve  the  equations  is  arranged  as  a  "cascade"  process.  This 
is  done  to  provide,  wherever  possible,  the  latest  values  for  quantities 
appearing  in  an  equation  to  be  integrated. 


This  technique,  when  applied  to  a  second  order  linear  oscillator  of  the 
form 


X 


=  O 


has  the  stability  properties  of  a  centered  difference  scheme.  A  rigorous 
stability  analysis  cannot  be  made  of  the  technique  applied  to  the  highly 
nonlinear  system  given  in  equations  (6-44)  to  (6-48).  It  1b,  however, 
hoped  that  the  stability  properties  demonstrated  on  the  linear  system  will 
be  retained.  Furthermore,  unlike  centered  difference  schemes,  all  possible 
up-to-date  information  is  used  in  the  solution  of  each  equation.  It  is 
therefore  hoped  that  the  problem  of  diverging  solutions  generated  by  alternate 
time  steps  mentioned  by  Haltirier  in  his  discussion  of  centered  difference 
schemes  (reference  4)  can  be  avoided. 


The  variables  are  obtained  in  the  following  order:  d  u~  u/  CT  / 

;  Y  >  )  ' 

Several  factors  influence  the  computational  sequence. 


. and  are  required  in  the  equations .  We  will 

thus  have  to  integrate  over  sigma  to  obtain  ^%^from  equation  (6-47)  and 


Note  that  both 


then  integrate  over  time  t^  obtain  ^  for  forming  Since  (f 

independent  of  6~,  we  solve  for  and  ^only  during  the  computations  for 
the  first <Tlayer. 


is 
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Also,  note  that  and  /jt  both  appear  in  equation  (6-48),  the 

equation  for  fT .  Equation  (6-47),  the  equation  for  ^%^als,  contains 
V<T'/'  •  In  order  to  have  consistent  ^-dependent  quantities  in  equation 
(6-46),  we  iterate  on  equation  (6-47).  First  equation  (6-47)  is  formed 
and  integrated  for  which  is  in  turn  integrated  for<^. 

is  computed  and  equation  (6-47)  is  re-evaluated  and  integrated  for 
The  latest  value,  of  and  the  value  of  ^%)^are  then  used  in  equation 

(6-48)  to  compute  Thus  the  value  of  V^used  in  the  equation  for  (p 

is  the  one  used  to  compute  the  value  of  which  also  appears  in  equation 

(6-  <  7) . 


In  equation  (6-45),  the  equation  for  l/cr  ,  there  is  a  triple  product, 
namely  (^)  (/ -  This  is  the  only  triple  product  appearing  in  the 

equations.  If  aliasing  errors  are  to  be  avoided  for  this  triple  product, 
a  much  finer  physical  grid  spacing  is  necessary  than  is  required  for  any 
of  the  other  terms  appearing  in  the  equations.  Rather  than  accept  the 
computational  penalties  imposed  by  the  finer  grid,  the  following  computational 
scheme  is  used.  The  product  t^.V^is  formed  in  physical  space,  transformed 
to  spectral  spacejind  back  to  physical  space  before  being  incorporated  in 
the  equation  for  This  assures,  given  proper  grid  spacing,  that 

truncation  of  the  spectral  expansion  is  done  to  avoid  aliasing. 


The  sequence  of  stepB  for  one  complete  integration  of  equations  (6-44) 
to  (6-48)  can  be  outlined  as  follows: 

Form  ^  nd  Integrate  for  ^  to  form 

Form  and  integrate  for  6”  .  Form  p 

Form  the  integral  equation  for  using  the  Cr-  and  c^.Pjust 

obtained  and  integrate  for 
Integrate  for  (f/  and  form 

Form  the  integral  equation  for  using  the  new  and 

Integrate  for 


o 

o 

o 

o 

o 


o 

o 


& 


.  and  integrate  for 
G* 


Form  the  equation  for 
Form  the  product  l r  perform  the  transforms  from 

physical  space  to  spectral  space  to  physical  space 
Form  ^(T^J^and  integrate  for  'T/<p- 
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IV  Computational  Considerations  of  the  Transform  Method 


This  section  deals  with  the  computational  requirements  imposed  by  the 
Transform  Method.  Most  of  these  requirements  are  imposed  to  avoid 
"aliasing"  in  transforming  from  physical  grids  to  spectral  coefficients. 

In  subsection  IV. 1  aliasing  is  defined.  General  solutions  are  outlined  in 
IV. 2,  and  particular  results  for  VSH  modeling  are  presented  in  subsection 
IV.  3. 


IV. 1  Aliasing 


Aliasing  arises  when  dealing  with  functions  that  are  expanded  in  series 
of  orthogonal  functions  such  as  Fourier  Series,  orthogonal  polynomials, 
scalar  or  vector  spherical  harmonics,  etc.  We  will  call  such  functions 
"Harmonics"  in  a  general  sense.  Suppose  a  function  F(x)  of  a  variable 
x  is  expanded  in  series  of  harmonics,  which  we  will  denote  .  We  assume 
the  are  orthogonal  over  some  interval  (0,L)  in  x;  then: 


(1-1) 


C 


n 


where 


r  l 

f  (x)  cpn  (x)  cl  X 


(1-2) 

(1-3) 


5 

In  order  to  obtain  the  values  of  f(x),  we  can  evaluate  equation  (1-1). 

Conversely,  given  all  the  values  of  f (x) ,  we  can  obtain  the  values  of  all 

the  coefficients  c  . 

n 

A  case  of  practical  interest  arises  when  f(x)  can  be  represented  in  a 
finite  number  of  expansion  coefficients,  i.e.,  when 


■  "1  Cn  f»  <*> 
n  ■  o 


if 


f(x) 


(1-4) 


In  this  case,  we  will  say  that  f(x)  is  band-limited.  It  is  well  known  that 

in  this  case,  to  evaluate  the  N  coefficients  we  do  not  in  general  require 

the  values  of  F(x)  at  all  > .  Instead,  given  a  certain  number  of  samples  of 

f(x)  (not  necessarily  equidistant)  we  can  obtain  the  coefficients  c  .  These 

n 

coefficients  will  reconstruct  the  function  f  (x)  exactly  (by  equation  1-1). 
The  number  of  samples  required  depends  upon  the  form  of  the  f(x). 

On  the  other  hand,  if  we  are  given  less  than  the  above  critical  number 
of  samples,  the  that  are  obtained  from  eg.  (1-2)  will  not  be  correct. 

The  "lower  frequency"  coefficients  (n=  0,  1,2...)  will  be  in  errov .  This 
condition  is  known  as  aliasing.  It  applies  to  all  spectral  modeling 
activities  since  they  all  work  with  truncated  series,  i.e.,  functions 
that  are  forced  to  be  band-limited. 

The  value  of  N  needed  to  represent  f(x)  in  1-4  will  be  called  the  Harmonic 
Content  of  f (x) . 

IV. 2  Avoiding  Aliasing  Effects 

In  order  to  avoid  aliasing  completely,  two  general  approaches  can  be 
followed: 

(1)  Ensure  that  enough  data  samples  are  in  hand  to  compute 
coefficients  accurately, 

(2)  Guarantee  that  all  functions  arising  in  modeling  activities 
are  properly  band-limited. 

Both  approaches  must  be  followed  in  practical  modeling.  In  general,  most 
problems  arise  in  nonlinear  operations.  For  example,  if  two  sinusoidal 
terms  of  a  given  frequency  are  multiplied  together,  the  result  is  a 
sinusoidal  term  of  double  the  frequency.  In  general,  nonlinear  operations 
raise  the  harmonic  content,  in  the  sense  that  more  terms  in  the  series  of 


*  '2 


type  (1-4)  are  necessary  to  express  the  result  exactly,  even  If  the 
original  terms  were  band-limited.  For  example,  a  division  introduces 
an  infinite  harmonic  content. 

Since  the  transform  method  consist’  of  transforming  variables  to  a 
physical  grid,  and  performing  the  nonlinear  operation  on  this  grid, 
the  general  outline  of  the  procedure  that  must  be  followed  in  spectral 
models  with  the  transform  method  is: 

(1)  Tailor  the  computational  procedure  such  that  nonlinearities 
which  yield  a  harmonic  content  higher  than  the  model's 
truncation  limits  are  minimized.  For  example,  avoid 
division,  reduce  number  of  products,  etc. 

(2)  Once  the  above  is  done,  evaluate  each  non-linear 
operation  to  determine  the  maximum  harmonic  content 
of  all  the  nonlinear  operations. 

(3)  Determine  the  size  and  spacing  of  the  physical  grid  such 
that  the  maximum  harmonic  content  introduces  no 
aliasing. 

When  the  procedure  is  performed,  we  must  know  the  original  harmonic 
content,  the  types  of  nonlinearity,  and  the  effect  of  each  nonlinearity 
on  the  harmonic  content.  For  example,  the  nonlinearities  introduced  by  the 
advection  term  (V.V)  v  are  known  from  its  analytic  formula,  Section  I. 

IV. 3  Avoiding  Aliasing  in  VSH  Modeling 

Two  general  type"  of  functions  appear  in  VSH  modeling:  polynomials  and 
circular  functions.  In  addition,  the  radial  functions  chosen  may  be  of 
another  type.  But  for  the  modeling  discussed  in  Section  III,  we  have  speci¬ 
fically: 


1)  Circular  functions  e.U" in  the  vector  and 
scalar  spherical  harmonics. 

2)  Associated  Legendre  polynomials  P  £  in  vector  and  scalar 

sperical  harmonics. 

3)  Legendre  polynomials  P°  -  P^  in  radial  functions. 

All  of  these  functions  are  orthogonal,  and  we  can  find  the  coefficients 
by  applying  (1-2).  It  is  well  known  that: 

1)  For  polynomial  type  functions,  the  best  method  of  reconstruction 
is  Gaussion  integration.  With  this  method,  a  polynomial  of 
degree  K  ^  2N-1  can  be  integrated  exactly  over  an  interval 
(-1,1)  with  N  sample  points  located  at  the  zeros  of  the 
Legendre  polynomial  of  degree  N. 

2)  Circular  functions  do  not  use  Gaussion  integration.  The 
K  coefficient,  K^(N-l)/2  can  be  recovered  with  N 
equally  spaced  points  over  the  interval  equal  to  the 
period  of  the  fundamental,  (0,2^). 

using  these  basic  relationships,  we  can  derive  the  rules  which  are  used 
to  compute  grid  sizes  for  the  transformation.  The  rules  are  given  for 
h  specific  function  type  below.  These  are  based  on  a  maximum 
'harmonic  generating'  mechanism  equivalent  to  the  product  of  two  variables. 

a.  Radial  functions  P,  (s),  N  «•  0,1,2,  ...  N.-l 

1* 

-1  <  €  <  1 

The  number  of  layers  required  is  given  by 


P3  (s)  ) 


For  example  If  N  *  4  (i.e.  we  use  P  (s) . 

then  (3.4-2)/2  -  5 

Points  have  a  Gaussian  spacing. 

b.  Colatitude  functions  P®  (C0Sj&  ) 

0,1,2,...  L-l 

9<  0  <  rP  ;  -1  <  COS  &  <  1 

The  number  of  ©  -  points  is  given  by 

N0  >  31^1  (3-2) 

? 


For  example  if  L-10,  N  >.  15 

e 

(The  model  uses  *  16  for  symmetry  reasons) 

Points  have  a  Gaussian  spacing. 

2.  A ACP 

c.  Longitude  functions  £ 

l  M.|  <  0,1,2,...  M-l 


The  number  of 
For  example. 


cf  -  points  is  given  by 
N<p  '>+)*■!-  L 
lf  M-10,  Xy  ^  39 


The  model  uses  40  for  symmetry  reasons.) 

The  points  must  be  equally  spaced  on  (0,  27?  ), 

If  a  Fast  Fourier  transform  were  used,  the  above  would  need  some 
modification. 


A  few  additional  comments  are  in  order.  In  the  model  equations  in  Section 
III,  there  are  two  divisions.  One  is  given  by  c/r  terms.  This  case  is 
acceptable,  because  when  (7  is  developed  in  Legendre  polynomials,  then 
0  is  a  factor  of  and  therefore  the  division  is  exact.  The  c-  _r 
division,  in  the  forcing  fields^ill  in  fact  introduce  some  aliasing. 

In  general  any  forcing  field  or  boundary  condition  causes  aliasing. 

If  this  aliasing  proves  empirically  troublesome,  then  the  grid  size  will 
have  to  be  expanded,  which  means  that  additional  coefficients  would  be 
used. 


One  triple  product  is  present  in  the  model  in  Section  III.  Aliasing  in 
this  product  is  avoided  by  using  the  product  computations  twice  and 
truncating  the  intermediate  result. 

IV. 4  Verification  of  Transform  Method 


A  set  of  utility  programs  for  manipulating  scalar  and  vector  spherical 
harmonics  have  been  programmed  in  the  APL  language.  These  include: 


Clebsh  -  Gordan  coefficient  computation 

Racah,3-J  and  6-J  symbol  computation 

Product  of  Scalar  harmonics  (Y?  Y^1,) 

Scalar  times  vector  harmonics  (YV,  t“  ) 

M  L  M,JL 

Vector  harmonic  dot\^product  (T^  . 

Vector  harmonic  cross  product  (t”  I,,, ,, 

JL  x  J  L  , 

Gradient  &  Laplacian  of  scalar  harmonics 

Divergence,  Curl,  Laplacian  of  vector  harmonics 

Advection  (T*1-.?  ,  Jl' 

JL  '  TT,L' 


') 


These  routines  are  coded  to  operr.te  entirely  on  harmonic  expansion 
coefficients. 


The  programming  system  computes  the  expressions  for  the  derivatives  of 
the  forecast  variables  in  physical  space.  The  derivatives  (physical)  are 
then  expanded  in  spectral  coefficients.  Each  term  in  the  expansion  of 
the  physical  derivative  is  the  derivative  of  the  corresponding  term  in 
the  expansion  of  the  physical  variable. 

To  test  the  code  that  does  the  spectral  to  physical  transformations  and 
to  verify  the  validity  of  the  transform  approach,  a  series  of  comparison 
runs  were  made. 

In  these  runs  the  scalar  and  vector  fields  in  the  model  were  initialized 
to  several  different  sets  of  coefficient  values.  The  program  system  was 
then  run  through  one  cycle  that  included  the  following  steps: 

a.  Transform  variable  fields  to  physical  space 

b.  Compute  spatial  derivatives  in  coefficient  space 

c.  Transform  results  of  (b)  to  physical  space 

d.  Form  a  comprehensive  set  of  non-linear  terms  in  physical  space 

e.  Transform  results  of  (d)  to  spectral  space. 

The  output  of  (e)  was  compared  to  the  results  of  APL  runs  where  the  various 
non-linearities  were  formed  using  the  fundamental  relationships  of  vector 
and  scalar  spherical  harmonics  on  a  term  by  term  basis  on  the  initial 
field  coefficient  values.  Both  the  program  system  and  APL  operate  in 
double  precision  (approx. 16  decimal  digits  resolution) .  In  all  cases 
there  was  agreement  to  the  thirteenth  digit. 


V. 


Radial  Functions  Analysis 


V.I  Introduction 

Expansion  of  vector  or  scalar  variables  In  vector  or  scalar  spherical 
harmonics  allows  us  to  separate  the  variables  of  a  partial  differential 
equation  In  the  sense  that  the  angular  part  can  be  Integrated  out,  leaving 
a  system  of  partial  differential  equations  In  the  harmonic  coefficients. 

If  the  problem  la  not  time-dependent,  the  equations  become  ordinary  dif¬ 
ferential  equations  In  r  but  In  the  case  of  a  time-dependent  problem, 
the  spectral  coefficients  are  functions  of  r  and  t.  By  expanding  the  co¬ 
efficients  again.  In  some  aet  of  suitable  functions  of  r,  the  problem 
can  be  reduced  to  a  aet  of  ordinary  differential  equations  in  time. 

The  approach  could  be  further  extended  to  expand  the  remaining  coefficients 
in  functions  of  time,  leaving  us  with  a  purely  algebraic  problem.  As  an 
example,  consider  the  partial  differential  equation 

dv  -  Vxv  ;  V  -  (1_ 

<3t 


In  order  to  solve  this  equation,  we  begin  by  expanding  the  vector  field 
v  in  vector  spherical  harmonics: 

h,L 

using  the  relations  for  curl  gives  (1-2) 
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t  ivr ,l 
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S-  i 


(1-3) 


and  equating  tha  coefficient*  gives  the  set  of  partial  differential  equations 
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,  bL  9  and  are  functions  of  r  and  t. 

In  this  particular 

case,  the  exact  functional  fora  can  be  determined  by  formal  separation  of 

variables  In  equations  (1-4)  through  (1-6).  But  suppose,  as  la  the  case 

»  with  more  complicated  problems,  that  we  are  unaware  of  the  exact  functional 

fora.  Then  we  can  Illustrate  the  expansion  process  by  choosing  a  set  of 

MM  M 

radial  functions  and  expanding  the  a^  ,  ,  and  coefficients.  Suppose 

we  choose  Fourier  Series  as  basis  functions  (because  of  their  simplicity) . 
Let  us  therefore,  expand  In  a  truncated  series  of '2N  Fourier  terms 


(nn 
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iter- 
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That  Is,  the  new  unknown  coefficients  are  functions  of  time  alone.  The  co¬ 
efficient's  radial  dependence  Is  now  explicit,  so 


(?,  t)  _  y  cb?k  ( t )  eztln 
dt  K  c,t 


(1-10) 
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with  similar  relations  for  bJJ  and  .  Then  using  (1-7)  -  (1-10)  in  (1-4) 


(1-6)  we  find 
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as  the 

e  functions 

are  orthogonal. 

K, 

L/C 


(1-14) 

(1-15) 
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For  any  given  rt  equations  (1-14)  -  1-16)  ere  now  ordinary  differential 
equations  in  tine;  each  value  of  M  and  L  yields  N  equations,  each  of  the 
type  (1-14)  -  (1-16).  Thia  set  happens  to  be  uncoupled  because  the  example 
chosen  is  very  slnple. 


lltr 

A  better  choice  of  radial  functions  in  this  case  is  the  set  of  functions  e  /r 
since  they  happen  to  be  eigenfunctions  of  the  radial  operation  ,  with 
eigenvalue  (lk).  In  this  case  we  would  have  obtained 
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(1-18) 


(1-19) 


In  fact,  the  appropriate  radial  functions  in  this  case  are  spherical  Bessel 
functions  (1).  This  can  be  ascertained  either  by  substitution  or  by  formal 
separation  of  variables  (1-4)  through  (1-6).  However,  the  intent  here  is 
to  Illustrate  the  procedure  and  not  to  particularize  results. 

Although  the  Fourier  Series  chosen  do  not  lead  to  as  simple  a  form  as  the 
functions  eikr/r,  they  have  the  advantage  of  behaving  well  under  product- 
type  nonllnearltles,  since  products  can  be  replaced  by  sums. 

The  preceding  discussion  Illustrates  some  of  the  desirable  properties  of 
radial  function  sets.  These  properties  will  be  discussed  In  more  detail  In 
the  following  subsection.  Subsection  V.3  will  review  some  candidates  in 
more  detail.  In  the  light  of  the  desirable  properties  discussed  In  V.2. 
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V.2  Desirable  Properties  of  Radial  Functions 


V.2.1  Orthogonality 

It  la  flrat  of  all  desirable  that  the  ^cabers  of  the  set  of  radial  functions 
chosen  (henceforth  "candidates")  be  mutually  orthogonal.  This  Is  obviously 
advantageous  froa  the  analytic  point  of  view.  From  the  computational  point 
of  view,  orthogonality  la  desirable  because  It  peralts  the  calculation  of 
expansion  coefficients  without  "aliasing"  effects  (see  Section  IV).  Also,  the 
computation  of  expansion  coefficients  is  usually  simplified  by  the  orthogonality 
condition. 

It  Is  also  desirable  that  the  functions  all  be  orthogonal  over  the  same  Interval 
and,  for  Mteorologlca^  purposes,  that  the  interval  up  conveniently  Into  the 
radial  dimension  of  the  atmosphere  (as  a  counter-example,  Bessel  functions  are 
orthogonal  only  over  zeros  cf  adjoining  functions). 

V.2. 2  Derivative  Properties 

The  radial  differential  operator  which  always  occurs  In  VSH  expansions  is 


L 


&  j  1 ,  2t- 


•  a 


(2-1) 


Thus,  it  is  desirable  that  the  candidates  have  the  closure  property  under  this 
operator,  i.e.,  that  the  operator  applied  to  a  member  of  the  set  yield  a  linear 
combination  of  members  of  the  set.  That  is,  if  <fi  is  a  member  of  the  candidate 
set, 


4  <4  (r>  -  1  c‘< 

where  the  limits  of  the  summation  depend  on  the  set  of  functions  and  L.  The  Ck 
are  coefficient  which  also  depend  on  the  set  of  functions  and  L.  This  property 
permits  us  to  analyze  the  "harmonic  content"  of  the  derivative  expressions  and 
is  therefore  useful  computationally  in  estimating  aliasing.  As  an  example,  if 
is  a  polynomial  of  degree  K  in  r,  its  derivative  is  a  polynomial  of  degree 
K-l.  It  csn  obviously  be  expressed  as  sum  of  at  most  K-l  <^pn’s,  nfK-1. 

If  we  only  require  one  term  in  (2-2)  then  ^i  is  an  eigenfunction  of  DL  .  Usually 
it  is  not  possible  to  find  simultaneous  eigenfunctions  of  all  the  's  appearing 
in  a  given  problem,  so  we  must  approximate  by  using  a  truncated  expansion. 

If  the  candidate  is  not  closed  under  the  D  operator,  then  at  least  we 
should  require  closure  under  the  d/dr  operation  to  estimate  aliasing  effects. 


(2-2) 


V. 2 . 3  Product  Rule 

In  many  applications  we  must  deal  with  product  terms,  e.g.,  (V.ty  V.  One 
of  the  advantages  of  VSH  ii  that  products  of  VSH  or  SSH  are  reduced  to  sums. 
This  property  is  desirable  both  analytically  (because  it  simplifies  the 
equations  as  a  rule)  and  computationally,  because  it  is  possible  to  consider 


N 


V 


aliasing  ef facts.  What  we  require.  If  not  an  explicit  product  rule,  la  at 
least  the  knowledge  of  the  harmonic  content  of  the  product  of  candidate 
■embers.  Thus,  for  example,  In  a  Fourier  elne  series  the  product  of  any 
two  terms  Is 

Afh  ±in  tn  x  ■  Sn  s»> a  *  =  (  ~  &  (>»+”)>  +  V  ) 

(2-3) 


Thus,  the  harmonic  content  of  the  product  of  two  sine  series  Is  the  sum  of  the 
harmonic  content  of  each  series.  Similarly,  the  product  of  two  polynomials 
of  degree  K  Is  a  polynomial  of  degree  2K. 

V.2.4  Applicability 

This  criterion  Is  difficult  to  define.  Some  function  sets  which  meet  all  the 
criteria  -  orthogonality,  differentiability  closure,  product  rule,  etc.,  are 
not  suitable  candidate  sets  for  reasons  not  Immediately  clear  in  the  light  of 
previous  discussions.  For  example,  the  fact  that  Fourier  series  are  periodic 
causes  problems  when  these  functions  are  used.  One  could  add  a  criterion  of 
non* periodicity  but  this  seems  too  sweeping.  The  spherical  Bessel  functions 
are  more  applicable  than  Fourier  Series  even  though  the  Bessel  Functions  do 
not  meet  all  the  nbove  criteria.  (They  lack  a  product  rule.) 

llie  main  criterion  In  deciding  applicability  Is  practical  experience,  coupled 
with  analysis  based  on  analytic  work  on  linear  problems  (see  previous  subsection) , 
or  oc  fits  to  real  data,  or  both.  A  third  consideration  Is  programming  difficulty 
and  evaluatlcn  time,  although  in  the  current  system  most  of  the  data  needed  for 


fitting  is  precomputed,  and  the  actual  fitting  reduced  to  vector-matrix 
multiplications.  Since  nonlinear  calculations  are  carried  out  at  grid  points 
using  the  transform  method,  the  actual  mathematical  formulation  of  the  product 
rules,  D_  operator  rules,  etc.,  is  not  as  Important  as  the  knowledge  of 
the  harmonic  content  of  the  formulas,  from  a  computational  point  of  view. 


V.3  Results  of  Radial  Function  Studies 

V.3.1  Summary 

A  number  of  function  sets  has  been  investigated  to  date.  A  complete  statement 
of  the  results  would  be  extremely  tedious,  since  it  would  involve  statements 
of  the  pertinent  properties  of  the  functions  in  question.  In  all  cases,  the 
functions  are  well  known  and  documented  elsewhere,  notably  references  (1), 

(2),  (3).  The  policy  adopted,  therefore,  is  to  provide  a  summary  statement  of 
the  virtues  (or  otherwise)  of  the  function  sets  in  the  light  of  sub-section 
V.2,  and  to  comment  briefly  on  the  summary  if  necessary.  Table  V.l  provides 
the  summary  and  the  subsequent  discussion  provides  the  commentary. 


Functions 

Orthogonal 

Derivative 

Closure 

Product 

Rule 

Applicability 

Powers  of  r 

No 

Yes 

Yes 

Diverge  at  r  -  °° 

Powers  of  1/r 

No 

See  (j.) 

Yes 

Bad  choice  based  on 
fits  to  data 

Fourier  Series 

Yes 

Yes  (2) 

Yes 

Problems  with 
periodicity  at 
end-points 

Legendre  Polynomials 

Yes 

Yes 

Yes 

Current  candidates 
problems  with  mapping 
interval  of 
orthogonality 

Spherical  Bessel 
Functions 

Yes  (3) 

Yes 

No 

Best  choice  except 
for  lack  of  product 

Associated  Laguerre 
Functions 

Yes 

Currently  being 
investigated. 

Table  V.l 


Radial  Function  Summary 


NOTES: 


(1)  If  one  truncate*,  then  tens  on  the  derivatives  of  Inverse  powers  of 
r  will  be  lost. 

(2)  Closed  under  ;  the  operator  Opr  +  1/r)  introduces  son  problems. 

(3)  Orthogonality  Interval  is  very  cumbersome. 

(4)  Actually,  radial  Eigenfunctions  of  the  one-electron  atom. 

V.3.1  Powers  of  r,  and  Inverse  Powers  of  r  or  1/r 

of  r,  l.e. ,  a  polynomial  In  r,  have  the  advantage  of  being  easy  to 
compute.  However,  they  are  divergent  for  large  r.  Pits  to  da^s  indicate 
that  these  functions  are  not  particularly  applicable.  Static  linear  problems 
sometime*  yield  powers  of  r  as  radial  eigenfunctions.  Por  example,  the  solution 
of  the  equation 

F'  Cp  (r,  9,<P)  =■  O  (3.!) 


(Laplace's  equation)  can  be  written  as 
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with  A(  and  B,  constants.  However,  no  dynamic  problem  (i.e.,  one  Involving 
did  t  of  any  variable)  this  far  attempted  has  yielded  powers  of  r  or  inverse 
powers  of  r  as  a  solution. 


V.3.2  Fourier  Series 

Fourier  series  hav*  been  used  to  fit  data  with  variable  results.  In  some  cases 
good  fits  ar-  obtains,  but  in  others  there  are  problems.  These  problems  arise 
from  the  fact  that  the  Fourier  functions  are  periodic  but  the  atmosphere  is  not. 
Any  fitting  process  must  make  some  assumptions  as  to  how  the  data  will  be 
repeated  outride  the  Interval  of  the  atmosphere. 

As  an  exarole,  one  ccn  assume  the  daM  repeat  symmetrically  with  odd  functlc  a 
predominating,  in  which  case  a  sine  series  is  obtained.  In  each  case,  however, 
the  data  must  be  examined  to  determine  the  proper  assumptions,  or  ei£«  the 
series  may  not  fit  well  over  the  prime  interval.  The  resulting  complications 
have  led  us  to  discard  Fourier  series  as  radial  functions,  at  least  for  the 
present.  Otherwise,  Fourier  Series  are  a  good  choice  Insofar  as  the  criteria, 
of  the  previous  section  are  concerned. 


V.3.3  Legendre  Polynomials 

These  are  the  choice  of  at  least  one  other  spectral  modeling  group  (4)  (5). 

As  is  well  known,  the  Legendre  polynomials  are  an  orthogonal  set  of  poly¬ 
nomials,  with  derivative  closure  and  well-defined  product  rule.  Computationally 
they  are  very  well  behaved. 


The  main  drawbacks  of  Tegendre  Polynomials  are  that  (1)  they  are  unsuitable 
for  variables  that  tend  to  zero  at  the  upper  atmosphere,  such  as  pressure;  and 
(2)  they  do  not  arise  in  any  linearized  hydrodynamic  problem  considered  to  date, 
and  therefore  are  probably  not  "natural"  In  soae  sense. 

In  spite  of  these  drawbacks  they  are  being  used  as  radial  functions  In  the 
model  described  in  section  III. 5. 


V.Z  4  Spherical  Bessel  Functions 

Spherical  Bessel  Functions  are  In  many  respects  the  most  "natural"  choice 
for  radial  functions.  They  are  in  fact  Bessel  functions  of  half-integer 
order,  and  arise  in  the  solution  of  many  linear  problems.  For  example,  the 
sound  save,  equations  (section  II. 4),  equation  (1-1)  in  this  section,  and 
Schrodlnger'a  equation  In  quantum  mechanics  for  the  Eigenfunctions  of  a  particle 
In  a  "spherical  square  well"  potential  all  lead  to  Spherical  Bessel  functions 
for  the  radial  dependence. 


With  ref  »rd  to  the  D  operator,  it  can  be  shorn  from  the  recursion  relations 
in  ref.  (1)  that,  if  f„(r)  is  a  spherical  Bessel  function  of  order  n,  then 
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The  simplicity  of  this  result  does  not  extend  over  into  products.  Thus  far. 
It  has  not  been  possible  to  derive  a  product  rule,  much  less  to  show  that  the 
set  is  closed  under  products.  Products  appear  to  (a)  double  the  argument  of 
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the  function,  and  (b)  introduce  a  polynomial  in  l/r  which  is  not  a  Spherical 
Bessel  function. 

With  regard  to  orthogonality,  the  Bessel  functions  are  orthogonal  over  zeros 
of  adjacent  functions,  which  introduces  a  computational  problem. 

The  literature  on  these  functions  is  not  extensive,  which  has  also  delayed 
their  evaluation. 

For  all  these  reasons,  Spherical  Bessel  functions  remain  a  potential  choice 
but  have  not  been  used  to  date  in  implementations. 


V.3.5  Associated  Laguerre  Functions 

To  be  exact,  we  refer  to  the  radial  eigenfunctions  of  the  one-electron 

M 

atom  (6).  These  arise  in  the  solution  of  Schrodinger's  equations  for  a 
particle  in  a  l/r  potential  field.  These  functions  are  currently  under 
Investigation. 
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VI.  Data  Handling 

This  section  is  a  report  on  numerical  work,  fitting  publlrhed  data 
with  vector  spherical  harmonics. 

The  results  of  the  fits  provide  information  on  the  number  of  terms  needed 
to  provide  accuracy,  on  the  radial  behavior  of  the  atmospheric  fields, 
and  on  computational  problems  that  will  be  encountered.  The  principal 
results  encountered  are  summarized  below. 

o  Number  of  Terms.  In  terms  of  the  index  L  in  the  scalar  or 
vector  spherical  harmonics,  it  appears  ‘hat  L*10  to  L-20  is 
sufficient  to  fit  the  published  data,  with  ample  accuracy. 

Since  most  published  data  are  zonally  averaged,  it  is  diffi¬ 
cult  to  make  similar  claims  for  the  index  M.  (The  index  L 
determines  the  maximum  frequency  in  &  ,  the  index  M  the 
maximum  frequency  in  (p  ). 

o  Radial  Behavior.  No  overall  conclusion  to  be  drawn.  In 
general,  there  appear  to  be  periodicities  combined  with 
overall  trends  (e.g.,  pressure  approaches  zero).  The 
radial  oehavior  can  be  duplicated  with,  for  example, 

Legendre  Polynomials  or  Fourier  Series  (see  Section  V). 
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o  Computational  Problems.  No  major  computational  problems 
arise  on  the  fitting  of  data,  except  those  connected  with 
aliasing  (see  section  IV).  In  the  case  of  using  Fourier 
Series  for  radial  representation,  additional  problems  are 
encountered  (section  V).  All  computations  presented  herein 
were  performed  in  APL  vichout  any  difficulty.  When  fitting 
data  it  is  advantageous  to  use  a  grid  spacing  of  "gaussian" 
points  to  perform  the  necessary  quadratures.  Thus,  it  is 
necessary  to  reduce  the  data  to  such  a  grid  prior  to  fitting. 
In  practice  this  could  present  some  difficulties.  As  with 
all  data  reduction  schemes,  data  gaps  pose  a  problem. 

The  remainder  of  this  section  consists  of  tabulated  results. 
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ZOHALLY  AVERAGED  MEAN  GEOPOTENTIAL  HEIGHT 
1000  to  700  MB 
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-103.76 
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2.88 
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236.38 
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-70.57 
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191.81 
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152.82 
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178.70 

-53.38 

157.15 

61.39 

174.33 

-47.81 

131.80 

49.30 

148.70 

-36.48 

111.79 

43.28 

122.26 

-30.54 

83.73 

34.32 

94.67 

-21.54 

61.35 

25.93 

69.84 

-14.40 

39.02 

15.68 

42.93 

-5.06 

25.08 

10.79 

25.96 

-1.79 

11.19 

5  88 

12.39 

-0.50 

5.26 

2.80 

4.95 

10770.57 

-422.05 


191.66 

103.86 

101.70 


Program  System  Development 


VII. 1  Work  to  Date 

Program  modules  have  been  developed  to  perform  all  the  functions  necessary 
to  apply  VSH  to  partial  differential  equations  expressed  in  npherical 
coordinates.  These  Include  modules  to: 

o  Numerically  fit  spectral  data  to  obtain  estimates  of  radial 
dependence 

o  Form  radial  differential  operators  in  spectral  space  for  use 
in  forming  three-dimensional  3patial  differential  operators 
o  Form  three-dimensional  spatial  differential  operators  -  curl, 
divergence,  etc.  —  in  spectral  space 
o  Form  horizontal  spatial  differential  operators 

o  Transform  between  spectral  and  physical  sp; ce 
o  Perform  time  and  space  integrations  in  Bpectral  space 
o  Perform  numerical  quadrature  in  spectral  space 
o  Generate  tabular  data  for  the  transformation,  fitting  and 
quadrature  modules 

o  Generate  the  job  control  language  procedure  to  execute  the 
programs  and  for  describing  data  sets  to  be  used  by  them 

In  addition,  data  organization,  both  for  core  and  peripheral  storage,  has 
been  developed  to  support  and  inform  these  program  modules.  Finally, 
programs  to  control  the  sequence  of  execution  of  the  program  modules  and 
to  form  the  model  described  in  Section  3.7  have  been  implemented. 

Programs  have  been  developed  in  terms  of  the  general  system  described  in 
the  previous  report.  Capabilities  not  described  in  that  design  description, 
but  required  for  the  Mintz-Arakawa  implementation,  such  as  spatial  integration 
and  numerical  quadurature,  have  also  been  implemented.  In  all,  some  9,000 
lines  of  code  have  been  generated  since  January  of  1973. 

The  equations  described  in  Section  3.6  and  implemented  as  in  Section  3.7 
have  been  programmed.  Debug  runs  ire  being  made  on  an  IBM  System/370  Model 
155  computer. 
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VII. 2 


Lessons  Learned 


It  is  somewhat  iror.ic  that  a  program  system  that  was  designed  with  the 
intention  of  anticipating  change  should  find  it  difficult  to  accommodate 
certain  functions  not  clearly  foreseen  at  the  time  of  the  design.  But  work 
on  the  Mintz-Arakawa  implementation  has  poi  ted  up  several  areas  where 
rethinking  and  reworking  are  in  order. 

The  problem  of  aliasing  when  dealing  with  transformations  between  spectral 
and  physical  domains  was  not  fully  understood  at  the  time  of  the  original 
design.  Truncations  of  spectral  expansions  over  theta  and  phi  were  done 
to  avoid  aliasing,  but  the  need  for  a  consistent  radial  truncation  was  not 
appreciated. 

It  was  the  original  intent  of  the  system  design  to  isolate  the  model  in 
physical  space.  The  equations  would  be  formed  and  boundary  conditions 
applied  in  one  compact  and  related  set  of  program  modules.  There  are, 
however,  instances  where  portions  of  the  model  might  more  efficiently  be 
formed  in  spectral  space.  Certain  boundary  conditions  are  more  readily 
cast  and  appl  d  in  the  spectral  domain.  So  rather  than  isolate  a  given 
model,  it  may  be  more  advantageous,  even  from  a  research  view,  to  distribute 
its  elements  through  the  system. 

The  system  was  originally  conceived  for  solving  temporal  differential 
equations.  The  final  cast  of  our  Mintz-Arakawa  model,  involving  differential 
equations  in  ^  and  <T  has  made  it  obvious  that  spatial  differential  equations 
must  be  accommodated.  Also,  while  specific  integration  algorithms  were  not 
considered  during  the  design,  the  need  for  sequential  rather  than  simul¬ 
taneous  integration  schemes  has  made  the  control  flow  for  integration  more 
cumbersome  and  less  straightforward  than  originally  thought. 
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All  of  these  considerations  can  be  fit  into  the  current  system.  With  the 
exception  of  forming  model  equations  in  spectral  space,  they  all  have  been 
for  the  Miotz-Arakawa  implementation.  These  changes  have  not,  however, 
been  made  with  the  general ' :y  and  coherence  that  the  original  design  strove 
for.  The  key  to  the  latter  is  anticipation  of  functional  requirements  and 
careful  planning  fot  t  leir  needs.  This  is  particularly  true  for  data 
organization.  The  heart  of  the  current  system  is  the  organization  and 
content  of  its  data  storage  areas.  Many  of  the  above  functions,  while 
implemented  in  the  system,  use  data  areas  originally  designed  for  other 
functions.  This  is  clumsy  and  at  times  restrictive.  If  the  system  is  to 
accommodate  tht  newly  appreciated  range  of  capabilities  with  the  originally 
intended  generality,  ease  of  use  and  case  of  change,  a  design  review  is 
needed.  This  is  particularly  true  in  the  areas  of  data  needs  and  organization, 
and  high  level  functional  control  flow. 


